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ABSTRACT 


A  simulation  model  has  been  developed  to  provide  a 
broad-based  capability  for  analyzing  water  pollution  in  a 
complex  river  system.  In  principle,  one  powerful  management 
tool  is  systems  analysis,  wherein  mathematical  optimizing 
techniques  are  employed  to  effect  rational  tradeoffs  between 
several  competing  river  pollution  control  designs.  This 
objective  has  been  accomplished  by  integrated  utilization  of 
computer  based  simulation  and  dynamic  programming 
optimization  technique. 

A  discrete,  deterministic,  dynamic  programming 
algorithm  developed  in  this  work  is  shown  to  be  of  much 
use  in  the  areas  of  water  pollution  control  and  management. 
Given  a  set  of  conditions,  the  optimization  model  determines 
several  combinations  of  the  quantity  of  pollutents  and  the 
the  quantity  of  artificial  aeration  required  to  meet  a 
pre-specif ied  water  quality  standards.  The  objective 
function  is  the  minimization  of  the  sum  of  the  squares  of 
the  aeration  costs  and  the  costs  incurred  by  damaging  the 
quality  of  river  water  or  unnecessarily  improving  the 
system.  The  original  constrained  allocation  problem  is 
simplified  by  converting  it  to  an  unconstrained  one  via  the 
use  of  Lagrange  multiplier. 
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The  simulation  model,  including  the  dynamic  programming 
optimization  model  was  applied  to  the  North  Saskatchewan 
River.  A  typical  optimal  aeration  capacity  allocation 
policy  and  its  corresponding  water  quality  profile  for  the 
North  Saskatchewan  river  is  presented.  The  relationship 
between  the  total  available  aeration  capacity  and  Lagrange 
multiplier  is  also  developed  treating  weighting  factor  as 
parameters. 
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CHAPTER  ONE 


INTRODUCTION 

1 .  I  GENERAL  DESCRIPTION 

Since  the  turn  of  the  century,  research  that  has  been 
directed  towards  the  field  of  environmental  engineering  has 
resulted  in  an  extensive  array  of  technology  that  can  be 
employed  to  minimize  the  impact  of  man  on  his  surrounding 
environment.  Only  fairly  recently,  however,  have  certain 
researchers  directed  their  efforts  toward  a  mathematical 
description  of  a  variety  of  systems  both  natural  and 
man-made,  that  are  significant  to  environmental  engineering. 
Dynamic  and  steady-state  mathematical  models  which  have 
resulted  from  the  various  investigations  have  been  developed 
to  the  extent  that  they  can  be  employed  to  delineate 
solutions  to  real-world  questions.  Armed  with  a  descriptive 
and  structural  model,  an  environmental  engineer  can  evaluate 
design  and  operational  alternatives  as  well  as  control 
strategies  through  a  series  of  computer  simulations,  thus 
placing  him  in  a  preferred  position  in  decision-making. 

While  models  in  and  of  themselves  are  not  a  panacea,  within 
certain  bounds  and  limitations,  system  simulations  can 
provide  valuable  assistance  to  the  environmental  engineer. 

Although  models  which  have  evolved  during  the  past 
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decades  have  been  extensively  employed  for  description  and 
simulation  of  laboratory  research  results,  they  have  only 
occasionally  been  employed  for  the  design,  operation,  and 
control  of  man-made  systems  or  for  the  planning  and 
management  of  natural  systems-  There  is  no  doubt  that  the 
pollution  control  efforts  which  are  now  under  way  have 
already  borne  some  fruit.  Due  to  the  pronouncements  of  the 
more  lurid  popular  articles  during  the  last  2  or  3  decades, 
the  condition  of  some  of  the  most  heavily  used  streams  have 
improved  markedly.  Nevertheless,  serious  problems  of  water 
quality  still  confronts  the  world.  Their  solution  will  lie 
not  only  in  the  field  of  engineering  but  also  reguire  the 
very  best  efforts  of  economists  and  computer  technologists. 

The  decisions  of  a  control  engineer  can  be  no  sounder 
than  the  basic  information  and  insights  concerning  water 
pollution.  For  a  decision  maker,  therefore,  ready  access  to 
the  current  knowledge  in  many  areas  involved  in  pollution 
control  decision  is  a  matter  of  extreme  importance. 

1 . 2  Systems  Analysis  and  Pollution  Control 

The  electronic  computer,  one  of  the  modern 
technological  developments,  has  profoundly  changed  the 
scientific  world.  Its  ability  to  handle  complicated 
information  has  dramatically  expanded  environmental 
engineering  boundaries. 
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The  combined  pressure  of  growing  population,  rapidly 
expanding  industries,  and  increased  demand  on  agricultural 
lands  are  creating  serious  threats  to  the  quality  of  the 
environment.  For  years,  wastes  other  than  those  directly 
affecting  public  health  caused  little  concern.  But  as  the 
way  of  life  became  more  complex,  the  pressure  for  pure  water 
and  clear  air  touched  practically  all  segments  of  economy. 
New  and  increasingly  stronger  legislation  has  been  enacted 
in  an  attempt  to  protect  water  and  air  resources. 

Effective  pollution  control  requires  the  coexistence  of 
at  least  three  factors: 

1.  high  level  of  knowledge  and  dedication  on  the  part  of 
environmental  engineers  and  researchers, 

2.  a  sincere  Government  commitment  to  clean  waters, 
manifested  in  a  workable  regulatory  scheme. 

3.  public  understanding  of  the  nature  of  the  pollution 
problem  and  support  for  the  control  efforts. 

However,  the  discussion  of  the  factors  mentioned  above 
are  beyond  the  scope  of  this  work.  It  is  hoped  that  the 
materials  contained  in  this  work  would  be  both  informative 
and  useful  to  water  pollution  authorities.  However,  this 
work  is  a  mere  scratch  on  the  surface  of  the  body  of 
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knowledge  that  will  be  required  to  cope  adequately  with  the 
present  and  emerging  problems  of  water  pollution  control. 

It  is  hoped  that  the  optimization  by  dynamic  programming 
approach  via  the  lagrange  multiplier  idea  will  stimulate 
others  to  undertake  research  in  this  most  vital  and 
challenging  field  -  river  pollution  control. 

* • 3  Literature  Review . 

During  the  past  two  decades,  a  significant  effort  of 
research  has  been  devoted  to  the  study  of  pollution  control 
strategies  and  the  analysis  of  waste  water  treatment 
facilities.  As  a  result,  many  different  analytical 
methodologies  have  been  developed  using  linear  and  nonlinear 
programming  techniques.  However,  in  establishing  the 
managerial  goals  for  a  water  quality  management  system,  the 
deterministic  nature  of  quality  standards  has  played 
profound  roles  in  modeling  processes. 

The  methods  of  operation  research,  systems  analysis  and 
modern  mathematical  analysis  have  been  applied  with 
increasing  intensity  in  the  fast  few  years.  In  many  cases, 
the  methods  admit  direct  application  in  water  resources 
management.  Although  the  number  of  people  in  water 
resources  field  having  an  ability  in  system  field  is 
relatively  small,  these  few  have  made  substantial  progress. 
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The  first  major  effort  to  apply  operation  research  and 
systems  analysis  in  water  resources  studies  was  by  the 
Harvard  Water  Program,  which  was  started  in  1956  and 
culminated  in  1962  with  the  publication  of  a  text  book-like 
report  of  the  research  objectives,  concepts,  methods  and 
techniques.  Fiering  (1961,  1964,  1967),  Hufschmidt  and 
Fiering  (1966),  Thomas  (1963)  and  Matalas  (1967)  are  among 
the  many  workers  of  the  group.  The  subjects  are  all 
concerned  with  operation  research,  systems  and  simulation 
techniques  in  the  solution  of  water  resources  design 
problems. 

Another  sprinkling  of  papers  came  from  Northwestern 
University  where  Charnes  influenced  the  use  of  digital 
simulation  which  includes  the  possibility  of  using  effluent 
storage  and  low-flow  augmentation  as  control  variable. 

Lynn,  Logan  and  Charnes  (1962);  Logan  (1962);  Lynn  (1964); 
Deininger  (  1965)  and  Heaney  (1  968)  are  examples  of  work  from 
this  source.  The  Lynn  influence  has  since  spread  to  Cornell 
University  where  such  work  as  Liebman  and  Lynn  (1966), 

Loucks  and  Lynn  (1966),  and  Loucks,  ReVelle  and  Lynn  (1967) 
have  resulted.  The  objective  of  all  their  simulation  study 
was  to  develop  a  model  to  represent  a  sewage  treatment 
system  in  which  various  combinations  of  ’treatment*  elements 
and  operating  disciplines  were  used  to  measure  the 
effectiveness  of  the  system  for  given  configurations. 


' 

!  ,  -  •  ’  •  4 

.  - 

■ 

. 


6 


In  U.  S.  Public  Health  Service  study  of  the  Delaware 
River  Basin,  extensive  use  of  these  methods  of  analysis  was 
reported  by  Thomann  (1967).  Grantham  and  Schaake  (1971) 
have  presented  the  use  of  simulation  as  a  tool  to  study  the 
surface  water  hydrology  of  a  stream.  Kolo  and  Haimes  (  1973) 
have  extensively  used  simulation  for  planning  and  expansion 
of  regional  water  resource  systems  for  Ohio  river. 
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1 . 4  Source  of  Pollution 

Pollution  of  streams  results  from  roan's  activities  at 
home,  in  factories,  and  in  form.  At  home  water  is  used  for 
many  purpose  most  of  which  involve  adding  solid  content  to 
water.  The  detergent  problem  resulting  from  domestic  water 
use  creates  many  difficulties.  Unfortunately,  it  is 
biologically  nondegradable  and  relatively  unaffected  by 
conventional  waste  treatment  methods.  In  wastes  it  produces 
foam  in  rivers  and  interferes  with  normal  water  and  waste 
treatment. 

Significant  pollution  also  results  from  industries  such 
as  chemical,  meat  packing,  dairy,  canning  and  tanning.  All 
industries  use  large  amounts  of  water  for  processing 
products.  The  used  water  contains  large  quantities  of 
decomposable  organic  matters  such  as  milk  waste,  offal  from 
animals,  decayed  fruits  and  vegetables.  In  general,  the 
polluting  power  of  such  industrial  wastes  will  exceed  the 
polluting  power  of  wastes  resulting  from  a  strictly  domestic 
population  by  ten  or  even  twenty  times. 

Technology  is  ever  advancing  and  new  products  are 
developed  and  old  ones  abandoned  before  the  pollution 
significance  of  some  of  them  can  even  be  determined.  Today 
seventy  five  percent  of  drug  store  prescriptions  involve 
drugs  unknown  five  years  ago.  In  many  chemical  industries. 
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the  majority  of  sales  involves  products  unknown  two  or  three 
years  ago.  Thus,  the  source  of  wastes  and  their  complexity 
can  be  expected  to  increase  directly  with  the  growth  of 
industry,  particularly  the  chemical  industry. 

1 . 5  Effect  of  Pollution. 

For  years,  health  educators  have  preached  against  water 
pollution,  citing  pollution  as  a  major  threat  to  public 
health,  esthetics,  and  economic  use  of  water  resources.  The 
measure,  not  only  of  the  "cleanness"  of  a  river,  but  also  of 
the  capacity  of  a  river  to  absorb  pollutional  loads  placed 
upon  it,  is  the  amount  of  dissolved  oxygen  (D.O)  which  the 
river  contains.  The  absence  of  dissolved  oxygen,  combined 
with  industrial  and  municipal  pollution  produces  a  septic 
condition  which  not  only  prohibits,  the  recreational  use  of 
water,  municipal  water  supply,  the  presence  of  fish  life  but 
also  creates  a  taste  and  odor  problem  for  downstream  users 
of  river. 

Pollution  can  also  significantly  increase  the  costs 
involved  in  subsequent  uses  made  of  polluted  water.  It  can 
increase  the  cost  of  water  treatment  for  municipal  and 
industrial  use. 

But  how  do  you  go  about  evaluating  how  much  pollution 
is  too  much?  How  far  back  should  we  roll  pollution?  How  do 
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we  evaluate  the  reduction  of  water  use  value  against  the 
cost  of  eliminating  the  pollution?  Certainly  we  do  not 
expect  to  return  all  our  streams  to  *wild  stream*  status. 
This  whole  area  requires  much  coordinated  interdisciplinary 
study  so  that  decisions  can  be  based  on  facts,  not  on 
desires  only. 
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1 . 6  Water  Quality  Criteria 

In  view  of  the  large  number  and  variety  of  uses  to 
which  water  is  put,  a  number  of  different  standards  of 
stream  quality  have  been  developed.  The  water  pollution 
control  agencies  serve  to  outline  in  detail  exactly  what 
constitutes  pollution  in  each  case.  In  establishing  a 
qualitative  standard  for  receiving  body  of  water,  the  water 
pollution  control  agency  must  consider  the  current  and 
future  water  needs  of  the  entire  region.  The  river  water 
may  be  used  for  a  source  of  water  supply,  for  recreation, 
and  for  receiving  wastes.  The  standard  of  permissible 
pollution,  thereafter,  must  be  established  on  a  basis  that 
will  bring  to  the  public  which  inhabits  the  river,  the 
greatest  benefit  for  the  least  cost. 

The  most  important  indices  representing  stream  water 
quality  are  biochemical  oxygen  demand  (BOD)  and  dissolved 
oxygen  (DO)  content.  The  BOD  concentration  denotes  the 
amount  of  oxygen  required  for  complete  biochemical 
decomposition  of  the  organic  wastes  contained  in  water. 
Hence,  the  BOD  decomposition  indirectly  represents  the  total 
amount  of  organic  waste  (pollutant)  in  water.  The  DO 
concentration  represents  the  water  oxygen  content.  For  a 
natural  stream,  i.e.,  a  stream  unpolluted  by  waste 
discharge,  the  DO  concentrat ion  is  at  a  level  of  about  80  to 


. 


. 

. 

.  . 


90  percent  saturation  ( 7  to  14  mg/11  depending  on  stream 
conditions)  and  the  BOD  concentration  is  below  10  mg/1. 
Onder  normal  conditions,  a  DO  concentration  of  5  mg/1  is 
sufficient  for  maintaining  a  normal  acguatic  life  (Krenkel 
and  Parker  1969)  and  to  preserve  a  satisfactory  water 
quality  for  multipurpose  use.  When  the  level  of  DO  drops 
significantly  below  the  standard  due  to  a  heavy  pollutant 
load  discharged  somewhere  upstream,  artificial  aeration 
augmentation  can  be  uesd  to  raise  the  level  of  DO  to  a 
specified  level. 

1.7  Preview 


In  Chapter  2,  the  motivation  for  the  application  of 
digital  computer  as  a  tool  in  river  pollution  control  is 
presented  in  detail  and  a  few  of  the  well  known  computer 
languages  for  simulation  are  reviewed. 

In  Chapter  3,  the  details  of  the  proposed  method 
(Dynamic  Programming  Technique)  is  presented  with  a  typical 
example  problem. 

In  Chapter  4,  a  mathematical  model  of  the  system  and  a 
discretized  dynamic  programming  algorithem  is  formulated  for 
finding  the  optimal  allocation  policy.  The  original 


1  Defined  in  Appendix-C. 
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constrained  allocation  problem  is  simplified  by  converting 
it  to  an  unconstrained  one  via  the  use  of  Lagrange 
multiplier . 

In  Chapter  5,  details  of  application  of  the  model  to 
North  Saskatchewan  river  is  explained  along  with  the 
analysis  of  the  result. 

Chapter  6,  includes  the  summary  of  the  results  and  some 
proposals  for  further  research  in  this  area. 


.  -  6 
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CHAPTER  TWO 


ROLE  OF  COMPUTER  IN  POLLUTION  CONTROL 

2.  I  Introduction 

— — — - — —  T— — — — 

T he  use  of  simulation  as  a  tool  for  planning, 
construction,  and  operation  of  complex  systems  has  increased 
materially  in  the  past  decade,  principally  due  to  the  advent 
of  digital  computers.  Their  ability  to  handle  complicated 
information  have  dramatically  expanded  engineering 
boundaries.  In  so  doing,  it  has  simultaneously  created 
numerous  opportunities  for  the  application  of  mathematical 
ideas  and  methods  to  the  solution  of  traditional  scientific 
problems  and  made  possible  the  exploration  of  research  areas 
in  engineering  and  science  either  previously  unattainable  or 
undreamt  of. 

In  recent  years,  there  is  an  increasing  interest  in  the 
application  of  computers  to  water  pollution  control.  The 
successful  implementation  of  automated  process  control  in 
petroleum,  chemical,  steel,  and  other  industries  has 
strongly  stimulated  the  investigation  of  a  similar  control 
in  river  pollution.  Environment  administrative  personnel 
have  begun  to  realize  that  computer  simulation  of  river 
system  would  give  good  answer  to: 
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(a)  achieving  the  increasingly  stringent  requirements 
mandated  by  regulatory  agencies; 

(b)  combatting  the  rising  cost  of  power  and  qualified 
labor ; 

(c)  optimizing  intgrated  river  pollution  control  systems; 

(d)  generating  the  myriad  information  for  operations, 
management,  and  environmental  regulatory  agencies. 

Computer  applications  to  the  monitoring  and  reporting 
of  plant  functions  are  becoming  common  and  are  receiving 
widespread  acceptance  as  a  standard  practice.  This  fact  is 
so  fundamental  to  the  whole  approach  that  both  industries 
and  universities  are,  in  fact,  revising  classical 
engineering  practices  to  accommodate  the  computer.  This 
capability  can  be  employed  most  effectively  in  river 
pollution  control.  In  simulation,  the  computer  provides  a 
basis  for  finding  the  most  economical  design  as  well  as 
several  feasible  alternatives. 

2.2  Necessity  of  Computer  Simulation 

Any  river  pollution  surveillance  program  has  to  handle 
a  great  deal  of  data.  Efficient  operation  demands  that  the 
data  be  well  organized  and  be  available  for  use  when  it  is 
required.  For  this  a  good  information  storage  and  retrievel 
system  of  a  quality  and  magnitude,  which  only  a  digital 
computer  can  handle  is  essential.  An  example  of  what  such  a 
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system  can  handle  is  the  management  of  fifty  years  of  data 
for  nearly  twenty  four  waste  disposal  stations  (about  a 
million  values)  situated  on  the  banks  of  the  North 
Saskatchewan  River  near  Edmonton,  Alberta.  So  the  role  of 
the  digital  computer  is  important;  apart  from  the  fact  that 
it  makes  mathematical  computations  with  unbelievable  speed 
and  accuracy,  it  is  able  to  store  and  recall  mathematical 
formulas  and  data  in  copious  quantities.  The  whole  problem, 
mathematical  relationships  and  data,  must  be  clearly  and 
explicitly  defined.  The  computer  is  able  to  simulate  in 
fast  time,  generating  many  years  of  operational  data  in  a 
few  minutes. 

In  economic  terms,  the  computer  application  often 
achieves  a  very  low  cost  per  calculation,  while  requiring  a 
relatively  high  initial  cost.  It  should  be  recognized, 
however,  that  once  a  computer  program  has  been  developed, 
the  initial  cost  for  reuse  in  subsequent  projects  is  usually 
negligible  or  marginal.  Cost  per  calculation  is  not  the 
only  factor  which  should  be  considered  when  evaluating  the 
advantages  and  disadvantages  of  computer  usage.  The 
accuracy  of  calculation  and  the  total  amount  of  information 
that  can  be  generated  must  also  be  considered.  The  computer 
can  often  utilize  a  more  exact  set  of  equations  and  it  can 
process  information  approximately  one  million  times  faster 
than  a  human  being.  In  addition,  once  the  program  has  been 
debugged,  the  rate  of  mistakes  by  computer  is  negligible. 
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The  indispensability  of  the  computer  in  simulation  and 
optimization  of  river  pollution  control  model  will  be  quite 
obvious  from  chapters  three,  four  and  five  of  this  thesis. 

The  information  obtained  from  simulation  becomes  a 
comprehensive  source  of  operational  data.  This  information 
can  fulfil  management’s  needs  for  investigating  operational 
trends,  studying  operation  costs,  and  developing  data 
correlations  for  pollution  control.  Once  in  operation,  the 
computer  reporting  system  comprises  the  data  base  for 
further  improvement  of  the  model.  Such  operational  data  are 
vital  in  evolving  models  and  validating  them.  Further,  the 
computer  reduces  the  time  reguired  for  technical  and 
clerical  record  keeping,  minimizing  human  errors,  and 
enhances  data  accessibility.  The  data  available  in  the 
memory  of  the  computer  can  be  statistically  analyzed  and 
furnished  to  engineers  for  design  purpose.  Without  computer 
assistance,  preparation  and  collection  of  such  information 
may  become  too  tedious. 

The  value  of  simulation  lies  in  its  flexibility.  A 
simulation  system  can  be  operated  under  all  variations  of 
normal  and  extreme  conditions  and  can  be  checked  out  and 
evaluated  entirely  by  the  computer.  Many  guestions  can  be 
asked  during  simulation  and  the  validity  of  the  answers  is 
an  indication  of  how  well  the  model  represents  the  true 


system  or  prototype. 
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2 . 3  Methods  of  Simulation  in  Pollution  Control 

There  are  two  general  types  of  computers  used  in 
simulation  namely: 

Analog  (continuous-variable)  and  Digital  (discrete-variable) 
computers,  and  both  types  are  employed  in  simulation 
studies.  A  simulation  technique  is  classified  as  analog  or 
digital,  depending  on  the  way  the  problem  is  formulated  and 
the  type  of  computer  needed  for  the  formulation  and 
solution. 

2.3.1  Analog  Simulation . 

An  analog  computer  is  one  in  which  computation  is 
performed  by  varying  the  state  of  some  physical  elements  in 
which  the  variables  are  continuous.  When  analog  simulation 
is  employed,  all  parts  of  the  system  must  be  simulated 
simultaneously.  The  parallel  nature  of  analog  simulation 
therefore  makes  it  closer  to  system  behaviour  than  does 
digital  simulation. 

Each  component  or  function  of  the  analog  simulation 
system  must  be  simulated  by  one  or  more  components  in  the 
computer.  This  nature  of  analog  simulation  always  requires 
the  addition  of  more  computer  equipment  as  the  system  being 
simulated  becomes  more  and  more  complex.  For  small  systems 
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this  requirment  is  not  serious,  but  for  the  study  of  large 
systems  the  addition  of  more  simulator  elements  can  become 
expensive.  Moreover,  in  many  cases  further  additions  of 
simulator  elements  may  not  be  feasible  because  there  is  a 
practical  limit  to  the  number  of  simulator  elements  that 
will  work  together  satisfactorily. 

Furthermore,  the  accuracy  of  the  analog  computer  is 
limited  to  the  accuracy  of  the  physical  components  involved. 
Thus,  the  use  of  analog  simulation  in  river  pollution 
studies  becomes  attractive  when  a  special  part  of  a  general 
system  is  to  be  scrutinized,  and  analog  technique  becomes 
less  attractive  as  the  problem  under  study  becomes  more  and 
more  generalized. 

2.3.2  Digital  Simulation. 

Digital  simulation  is  characterized  by  the  use  of  a 
digital  computer.  Whereas  the  analog  computer  must  handle 
all  elements  of  the  simulation  system  simultaneously  (in 
parallel) ,  the  digital  computer  handles  elements  of  the 
simulation  system  one  after  another  (in  series)  .  Hence,  an 
increase  in  the  system  complexity  results  in  an  increase  in 
the  time  required  for  computation.  Although  accuracy  may 
not  necessarily  be  an  important  factor  in  the  simulation,  it 
is  possible  in  digital  simulation  to  reach  any  degree  of 
accuracy  desired  by  carrying  out  more  computations. 
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Moreover,  explosive  development  of  ever  faster  and  less 
costly  digital  computers,  coupled  with  better  software  and 
integration  routines,  seem  to  indicate  the  demise  of 
analog/hybrid  simulation,  except  for  some  interface 
equipment  needed  in  partial  system  tests  (e,g.,  with 
autopilot  components,  man-in-the-loop) .  As  a  case  in  point, 
training-type  flight  simulators  have  been  all-digital  for 
years.  As  a  rule  of  thumb,  a  1973  16-bit  minicomputer  can 
comfortably  simulate  a  three  dimensional  aircraft  engine 
equations  three  times  faster  than  real  time 
(G.  A.  Korn,  1973)  if  we  use  fixed-point  computation  and 
omit  fast  sub-system  such  as  hydraulic  servos. 

2.4  Simulation  with  FORTRAN 

In  discrete-event  simulation  the  language  used  must 
enable  the  system  designer  to  represent  a  complex  system 
conveniently  and  implement  comfortably.  The  representations 
could  involve  such  requirements  as  solving  of  equations, 
preserving  interrelations,  representing  decision  logics  and 
physical  characteristics  of  the  system.  The  languages 
employed  vary  from  Assembler  languages  through  the  general 
purpose  programming  languages  (FORTRAN,  PL/ 1 ,  ALGOL)  to  the 
general  purpose  Simulation  languages  (GPSS,  SIMSCRIPT, 
DYNAMO) .  In  RIPSIM  (River  Pollution  Simulation  Model) ,  the 
general  purpose  programming  language  FORTRAN  IV  was  used  as 
the  tool  over  the  available  general  purpose  simulatiom 
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languages,  notably  GPSS,  for  the  following  reasons: 

(a)  GPSS/360  is  available  only  on  larger  machines  so  that 
the  use  of  such  a  language  would  limit  the  scope  of 
the  model  and  hence  its  use  as  a  tool  for  the 
environmental  engineers. 

(b)  The  speed  of  execution  in  GPSS  tends  to  decrease  with 
the  growth  of  model  size  and  complexity;  and  as  a 
tool,  it  is  important  that  execution  be  as  fast  as 
possible  to  make  the  use  of  the  model  economical. 

(c)  It  is  easier  to  expand  or  decrease  the  size  of  the 


model 

through 

re dimensioning 

of 

the 

appropriate 

FORTRAN 

arrays 

than  it  would 

be 

for 

a  novice 

programmer  to  use  the  REALLOCATE  feature  in  GPSS  to 
expand  a  model  programmed  in  GPSS. 

According  to  Reitman ( 19 74)  the  general  requirements  of 
a  language  in  the  area  of  simulation  may  be  reduced  to  four 
basic  characteristics  as  follows: 

!.  Short-term  results; 

2.  Ability  of  system  to  represent  the  real  world; 

3.  Long-term  results; 

4.  Effort  required. 


1. 


Short-term  results.  The  system  designer  must  have  a 
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good  programming  background  in  the  language  before  he 
can  use  it  for  simulation.  But  FORTRAN,  being  one  of 
the  most  common  languages  gives  less  difficulty  as  far 
as  background  in  programming  is  concerned. 

2.  Ability  of  system  to  represent  real  world .  Almost 

any  real-world  condition  could  be  represented  in 
FORTRAN;  the  effort  goes  up,  however  with  complexity 
in  a  nonlinear  relationship. 

(a)  Logical  situations  in  the  model  can  be  represented 
without  difficulty. 

(b)  Mathematical  capability  of  the  language  is 
excellent.  There  are  numerous  special-purpose 
techniques  for  data  smoothing,  and  other  forms  of 
data  manipulations  which  are  both  available  and 
accessible  for  simulation. 

(c)  Maximum  model  size  is  completely  under  the  control 
of  the  programmer.  He  can  make  trade-offs  between 
storage  hierarchy  and  speed  of  execution. 

3.  Long-term  results .  The  most  important  advantage  of 
FORTRAN  in  this  area  lies  in  its  generality. 

(a)  Documentation  is  under  the  control  of  the 
individual.  There  are  aids  in  the  form  of 
cross-reference  files  after  the  individual  has 
thoroughly  set  up  his  comments. 
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(b)  System  designers  other  than  the  original  model 
developer  can  easily  follow  the  logic  and  detail 
of  the  simulation  with  less  effort. 

(c)  Computers  of  different  manufacturers  can  use  the 
same  higher-order  language  program.  Usually  there 
is  some  requirement  for  rework;  but  in  terms  of 
the  total  effort  involved,  this  would  be 
considered  minor. 

4.  Effort  required .  Less  effort  is  required  to  carry  out 
the  simulation  in  FORTRAN  than  it  would  be  for  the 
rest  of  the  General-purpose  programming  and  Simulation 
languages.  Several  programmers  can  work  on  the 
simulation  model  in  parallel  if  the  conventions 
governing  the  transfer  of  data  between  subroutines  are 
well  planned  in  advance.  In  this  respect,  a 
simulation  program  written  in  FORTRAN,  in  modular 
fashion,  can  easily  be  altered  by  anyone  with  a 
working  knowledge  of  FORTRAN,  and  in  particular,  parts 
of  the  program  can  be  altered  without  affecting  the 
whole  program. 

The  FORTRAN  meets  all  of  the  above  requirements  fairly 
well  and  further,  its  use  in  this  study  is  also  obvious  from 
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The  table  shows  the  relative  comparison  (Colella,  1974) 
of  various  languages  for  the  implementation  of  discrete 
simulation.  The  basis  for  these  weights  are  subjective,  but 
at  least  it  is  an  attempt  to  relate  these  considerations. 

The  languages  mentioned  above  are  merely  representative  of 
the  vast  number  that  are  available. 

In  the  table,  each  language  is  given  a  numerical  weight 
with  regard  to  a  specific  consideration.  The  weights  are 
interpreted  from  0  to  5  as  follows. 

0  -  not  applicable. 

1  -  very  poor. 

2  -  poor. 

3  -  fair. 

4  -  good. 

5  -  very  good. 


. 
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TABLE  2^1  Comparison  of  Languages  for  Discrete  Simulation. 
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FORTRAN 


PL/  \ 


GPSS 


Simscript 
_ I 


Execution 

.sjDeed 


Memory 

Utilization 


Availability: 
Oser’s  Library 
Personnel 
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Program  Tran¬ 
sferability  to 
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Application 
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Accuracy  of 
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ntation 


Simulation 

Structural 
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Programming 
Diaqostic  Aid 


Program  Docum¬ 
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CHAPTER  THREE 


DYNAMIC  PROGRAMMING  APPROACH  FOR  OPTIMIZATION 


3. 1  Introduction 


DYNAMIC  PROGRAMMING,  as  an  optimizing  procedure,  has 
been  around  for  a  few  years  now  and  perhaps,  the  main  reason 
that  it  is  still  with  us  arises  from  Bellman’s  inexhaustible 
supply  of  examples,  both  realistic  and  imaginary,  which  may 
be  analyzed  using  this  approach.  A  more  detailed  and 
comprehensive  coverage  of  dynamic  programming  can  be  found 
in  Bellman’s  three  books  (1959,  1962,  1965)  where  his 
examples  and  exercises  cover  almost  every  form  of 
optimization  problem.  Dynamic  programming  is  a  mathematical 
technique  often  useful  for  making  a  sequence  of  interrelated 
decisions.  It  provides  a  systematic  procedure  for 
determining  the  combination  of  decisions  which 
maximizes/minimizes  overall  effectiveness.  This  method  is  a 
general  type  of  approach  to  problem  solving,  and  the 
particular  eg-uations  used  must  be  developed  to  fit  each 
individual  situation.  Therefore,  a  certain  degree  of 
ingenuity  and  insight  into  the  general  structure  of  dynamic 
programming  is  required  to  recognize  when  a  problem  can  be 
solved  by  dynamic  programming  procedures,  and  how  it  would 
be  done. 
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Batches,  Haimens,  and  Hall  (1S69)  and  Morin  and  Esogbue 
(1971)  studied  the  optimal  sequencing  and  scheduling  of 
water  supply  projects  via  dynamic  programming,  minimizing 
the  present  value  of  the  total  cost.  Nainis  and  Haimes 
(1975)  expanded  the  Butcher,  Haimes,  and  Hall  (1969)  model 
for  the  optimal  scheduling  of  construction  of  multiple 
project,  multiple  location  of  water  supply  systems.  Hall, 
Haimes,  and  Butcher  (1972)  considered  the  construction  of 
interim  projects  with  limited  life  times,  such  as 
desalination  plants  via  dynamic  programming. 

Scarato  (1969)  discussed  the  time-capacity  expansion  of 
water  supply  systems.  Using  dynamic  programming  a  policy  of 
expansion  is  determined  to  minimize  the  present  value  of  the 
total  cost.  Rachford,  Scarato,  and  Tcfcobanoglous  (1974) 
apply  the  above  method  to  wastewater  treatment  plants. 

Kaplan  and  Haimes  (1975)  have  employed  the  dynamic 
programming  to  determine  the  optimal  schedule  of  expansions 
of  wastewater  treatment  plant.  Some  of  the  shortcomings  of 
the  above  method  are  overcome  by  a  method  proposed  by 
Riordan  (1976),  where  a  dynamic  programming  is  used  to 
determine  the  optimal  policy  of  capacity  expansion, 
minimizing  the  present  value  of  the  total  cost. 

Most  of  the  above  sources  are  concerned  with  water 
supply  and  waste  water  treatment  projects  and  dynamic 
programming  was  employed  to  determine: 
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1.  the  optimal  schedule  of  expansion  at  each  plant. 

2.  the  optimal  operating  policy. 

3.  optimal  allocation  of  computer  registers. 

4.  optimal  control  of  waste  water  treatment. 

5.  optimal  design  of  concrete  beams. 

In  this  work  a  discretized  dynamic  programming  is  used 
to  find  an  optimal  quality  control  of  river  water.  The 
objective  function  is  to  be  minimized  subjected  to  the  given 
total  available  aeration  capacity.  It  works  best  when  the 
number  of  decisions  at  any  stage  is  not  too  large.  Under 
this  condition,  it  can  handle  a  very  large  number  of  stages 
with  comparatively  little  difficulty.  The  philosophy  of  a 
sequential  decision  making  in  the  capacity  allocation  to 
each  of  the  series  of  aerators  along  the  stream,  the 
nonlinearities  in  the  DO  sag  profile  and  recursive  cost 
function,  suggest  the  use  of  dynamic  programming. 

3. 2  Characteristics  of  Dynamic  Programming 

The  principal  characteristics  of  any  problem  which  is 
amenable  to  solution  by  dynamic  programming  are: 

I.  The  problem  can  be  divided  up  into  stages,  with  a 
policy  decision  reguired  at  each  stage. 

The  stream  segment  under  consideration  is  divided 

into  N  reaches.  The  system  is  considered  to  be 


. 

>  I  L  » 

. 


. 


. 


■ 


28 


made  up  of  a  number  of  connected  segments 
(reaches)  in  which  each  reach  has  constant 
properties.  Reach  nodes  are,  therefore, 

established  at  every  major  or  significant  change 
in  the  system.  Significant  changes  are  the 
discharge  point  of  a  waste  load,  a  tributary 
joining  the  system,  etc.  A  long  section  of  a 
stream  having  acceptable  uniform  properties  should 
not  be  divided  into  two  or  more  reaches  because  of 
its  length.  The  policy  decision  at  each  stage 
(reach)  is  the  amount  of  aeration  V(i)  needed  for 
that  particular  stage. 

2.  Each  stage  has  a  number  of  states  associated  with  it 
and  each  state  corresponds  to  a  possible  condition 
which  the  decision  variable  might  take  at  any 
particular  stage. 

The  state  associated  with  each  stage  is  the  state 
variable,  D(i),  discretized  into  M  values  with 
C(s)  as  its  maximum  and  zero  as  its  minimum  value. 

In  general,  the  states  are  the  various  possible 
conditions  in  which  the  system  might  find  itself 
at  that  stage  of  the  problem.  The  number  of 
states  may  be  either  finite,  as  in  this  case,  or 
infinite. 


3. 


The  effect  of  the  policy  decision  at  each  stage  is  to 
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transform  the  current  state  into  a  state  associated 
with  the  next  stage. 

The  dynamic  programming  problems  can  be 
interpreted  in  terms  of  network.  Each  node  would 
correspond  to  a  state.  The  network  would  consist 
of  columns  of  nodesr  with  each  column 
corresponding  to  a  stage,  so  that  flow  from  a  node 
can  go  only  to  a  node  in  the  next  column  to  the 
right.  The  number  assigned  as  the  "distance" 
between  connected  nodes  can  be  interpreted  as  the 
amount  of  aeration  needed  to  go  from  one  stage  to 
another.  Now,  the  objective  would  be  to  find  the 
route  through  the  network  which  minimizes  the 
amount  of  aeration,  and  in  turn  minimizes  the 
objective  function. 

4.  Given  the  current  state,  an  optimal  policy  for  the 
remaining  stages  is  independent  of  the  policy  adopted 
in  previous  stages. 

Given  the  state  in  which  one  is  currently  located, 
the  optimal  aeration  capacity  policy  from  this 
point  onward  is  independent  of  how  we  got  there, 
since  the  knowledge  of  the  current  state  of  the 
system  conveys  all  of  the  information  about  its 
previous  behaviour  necessary  for  determining  the 
optimal  policy  henceforth. 
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5-  The  solution  procedure  begins  by  finding  the  optimal 
policy  for  each  state  starting  from  the  last  stage. 

6.  A  recursive  relationship  is  available  which  identifies 
the  optimal  policy  for  each  state  with  n  stages 
remaining,  given  the  optimal  policy  for  each  state 
with  n-  I  stages  to  go. 

Therefore,  finding  the  optimal  policy  when 
starting  in  state  s  with  n  stages  to  go  required 
finding  the  minimization  value  of  D (n)  .  This 
policy  would  consist  of  using  the  value  of  D (n) 
and  thereafter  following  the  optimal  policy  when 
starting  in  state  D(n)  with  (n-  1)  stages  to  go. 

7.  Using  the  recursive  relationship,  the  solution 
procedure  moves  backward  stage  by  stage  -  each  time 
finding  the  optimal  policy  for  each  state  of  that 
stage  -  until  it  finds  the  optimal  policy  of  N-stage 
system. 

8.  The  N-stage  optimal  policy  is  found  by  passing  forward 
through  the  system. 
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3 . 3  Philosophy  of  Dynamic  Programming  Technique 

How  does  all  of  the  above  characteristics  fit  in  with 
the  river  pollution  problem? 

Consider  a  river  basin  containing  one  major  stream  and 
several  industries  and  municipalities  located  adjacent  to 
the  stream  as  shown  in  Figure  3. I.  The  waste  discharged  to 
the  stream  is  assumed  to  be  characterized  by  its  BOD.  To 
meet  specified  DO  standards  for  the  stream  as  efficiently  as 
possible,  decisions  have  to  be  made  regarding  levels  of 
treatment  at  each  of  the  N  outfalls.  It  should  be  apparent 
that  the  river  system  can  be  divided  into  a  set  of  N 
discrete  reaches  or  stages  with  a  decision  being  required  at 
each  as  shown  in  Figure  3.2.  It  should  be  noted  that 
numerous  conditions  other  than  waste  discharge  can  and 
should  constitute  stage  or  reach  doundaries, 
e. g. , significant  tributary  flows,  change  in  temperature, 
change  in  flow  conditions,  etc.  These  cause  no  difficulty, 
however,  but  will  not  be  considered  in  this  application. 

The  state  of  a  reach  or  stage  consists  of  the  BOD  and 
DO  entering  at  the  upstream  end  of  the  stage.  There  are 
obviously  an  infinite  number  of  possible  states  or  incoming 
BOD  and  DO  combinations  for  each  stage.  The  range  of 
incoming  BOD  and  DO  is  divided  into  discrete  values  as  shown 
in  Figure  3.3. 
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FIGURE  3.1:  Schematic  Illustrating  a  Typical  River  Basin 
and  Waste  Outfalls. 


FIGURE  3.2:  Schematic  Illustrating  a  Typical  River  Basin 
and  Waste  Outfalls  as  an  N-Stage  Process. 
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FIGURE  3.3:  Schematic  Illustrating  Discrete  States  and 
Decisions  for  a  3-Stage  Process. 
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TABLE  3.1:  Number  of  Possible  Policies  as  a  Function  of 

Number  of  Stages  and  Treatment  Increments. 
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For  example,  if  incoming  BOD  could  vary  from  0  to  20  mg/1  at 

stage  n  one  might  select  to  use  values  of  0,  2,  4,  . , 

20  mg/1.  If  incoming  DO  could  vary  from  0  to  9  mg/1  at 

stage  n  one  might  use  0,  1,  2,  . ,  9  mg/1.  This  would 

result  in  110  possible  discrete  states  at  stage  n.  One  may 
use  increments  as  small  as  desired.  It  is  assumed  that  the 
BOD  and  DO  entering  the  upstream  is  known. 

At  each  stage  a  decision  regarding  treatment2  must  be 
made.  This  can  vary  from  0  to  9  mg/1.  The  consequences  as 
regard  cost  and  quality  must  be  investigated  for  the  whole 
allowable  range.  As  with  the  input  BOD,  discrete  levels  are 

used,  e.g.,  0,  5,  10,  . .  100  mg/1  and  DO  can  vary  from 

0,  I,  2,...., 9  mg/1.  Large  or  smaller  increments  may  be 
used.  The  relationship  between  number  of  stages,  treatment 
increment,  and  number  of  possible  N-stage  policies  is  shown 
in  Table  3.1.  There  are  an  infinite  number  of  possible 
policies.  Simple  preliminary  screening  can  reduce  the 
number  of  possible  policies  greatly.  Even  if  the  reduction 
is  several  orders  of  magnitude,  there  would  be  too  many 
alternatives  to  investigate  without  a  dynamic  programming 
approach. 

Transfer  functions  are  reguired  to  determine  the  BOD 
and  DO  leaving  stage  n  (and  entering  stage  n-1)  as  a  result 


2  Amount  of  aeration  reguired  at  each  stage. 
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of  being  in  a  given  state  (incoming  BOD  and  DO)  and  the 
decision  regarding  treatment  as  shown  in  Figure  3.4.  One 
may  use  transfer  functions  as  complex  and  realistic  as 
desired.  For  this  application,  the  basic  Streeter-Phelps 
(1925)  equations  for  BOD  and  DO  are  used. 

-k  t 

BODt  =  BODa  e  (3.1) 


(3.2) 


BODa  ,  the  initial  BOD  at  the  top  of  the  stage,  is  an  average 
of  the  incoming  BOD  (BOD  )  and  the  BOD  discharged  (BCDdis) 
at  the  outfall: 


BODa  =  £(BODin  ,  BOCdl^ 


(3.3) 


Since  BOD.  is  part  of  the  state  or  input  conditions  to  the 
stage  and  BODdis  is  a  discrete  function  of  the  treatment 
level.  Then  one  can  rewrite  the  above  equation  as: 


(3.4) 


BCD  in  =  f  (state,  decision) 


Similarly,  the  DO  should  be  an  average  of  the  incoming 


and  the  DC  in  the  waste  (BO^is 


=  f<DOin'  DOdis> 


(3.5) 
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FIGURE  3.4:  Schematic  of  a  Typical  Stream  Reach  or  Stage 

Illustrating  Transfer  Functions  for  BOD  and  DO. 
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Again,  DO in  is  part  of  the  state.  It  is  assumed  that  there 

is  a  negligible  amount  of  10  (DO  )  in  the  Waste.  But  for 

dis 

sake  of  generality,  the  transfer  equation  can  be  rewritten 
as: 

D0a  =  f (state,  decision)  (3.6) 

The  initial  DC  deficit,  DA ,  is  then: 

Da  =  Cg  -  D0a  =  f  (state,  decision)  (3.7) 

Now  by  substituting  B0DA  and  D^  from  equations  3.3  and 
3.7  into  the  BCD  and  DO  transfer  functions,  equations  3.1 
and  3.2,  the  EOD  and  DO  at  any  point  within  a  stage  is  a 
function  of  the  state,  or  input  BOD  and  DO,  and  the 
treatment  decision. 

BCDt  =  f (state,  decision)  (3.8) 

DO  =  f(state,  decision)  (3.9) 

Ey  determinig  the  DO  at  the  Sag2  point,  one  may  determine, 
for  any  state,  whether  the  DC  standard  has  been  violated  as 
a  result  of  various  treatment  level  decisions.  If  the 
standard  is  violated,  then  any  policy  which  include  that 
state  and  that  decision  would  he  infeasible  and  need  not  be 
considered  further. 


2  Defined  in  Appendix-B 
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When  If  equals  the  tiroe-of-flow  in  the  stage,  the 
output  EOD  and  DO  are  given  fcy: 


BOD 

out 


-k  Tf 

BODa  e  1 
A 


(3.10) 


k  BOD  -k  Tf  -k  Tf  -k0Tf 

DOout  =  Cs - -  >  +  D  e  2  (3.11) 

kl  -  k2 

Each  of  these  two  output  components  could  be  components 
of  the  state  for  the  downstrem  reach.  It  should  be  clear 
now  that  the  transfer  functions  enable  one  to  determine  the 
next  state  from  the  present  state  and  the  decision. 

Given  a  particular  state  or  input  BOD  and  DO,  the 
optimal  policy  for  the  downstream  reaches  is  independent  of 
the  upstream  decisions  that  may  have  lead  to  this  state. 

The  solution  process  begins  at  the  downstream  stage 
(stage  I)  and  moves  upstream.  Given  all  possible  states  or 
combinations  of  BOD  and  DO  which  could  enter  the  top  of 
stage  I,  one  simply  determines  the  lowest,  and  most 
economical,  level  of  treatment  which  will  prevent  a 
violation  of  the  DO  standard.  The  treatment  level  and  its 
cost  are  stored  with  the  state  for  future  use. 


At  stage  2  all  possible  states  or  combinations  of  BOD 
and  DO  entering  the  top  of  stage  2  are  considered.  For  each 


- 
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of  these  states,  the  entire  range  of  treatment  at  stage  2 
must  be  investigated.  The  guantity  which  is  to  be  minimized 
is  the  sum  of  treatment  costs  at  stages  2  and  I.  Selection 
of  a  treatment  level  at  2  allows  one  to  know  the  local  cost, 
as  cost  is  a  function  of  treatment  level.  If  the  level  of 
treatment  being  considered  does  not  cause  a  DO  standard 
violation  in  stage  2,  then  the  output  BOD  and  DO  are 
determined.  Knowledge  of  these  values  enables  one  to  go  to 
the  appropriate  state  in  stage  1  and  retrieve  the  optimal 

1- stage  cost.  The  entire  range  of  treatment  levels  for  each 
state  of  stage  2  is  used,  and  the  level  which  results  in  the 
least  2-stage  cost  is  the  optimal  treatment  decision  for  the 
given  state  of  stage  2.  Other  possible  states  of  stage  2 
are  investigated  in  a  similar  manner.  For  stage  2  (and  all 
other  upstream  stages) ,  the  optimal  treatment  level,  the 
optimal  total  cost  up  through  that  stage,  and  the  output  EOD 
and  DO  are  stored  for  each  state. 

When  one  reaches  stage  3,  the  information  about  stage  I 
will  not  be  used  directly,  as  it  is  incorporated  in  the 

2—  stage  information.  In  general,  all  that  will  be  needed  at 
stage  n  will  be  the  information  previously  obtained  for  an 
(n- h) -stage  system. 

When  the  upstream  (Nth)  stage  is  reached,  there  will  be 
only  one  possible  state  as  this  is  an  initial-value  problem. 
For  this  single  state,  all  treatment  decisions  will  be 
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considered.  The  quantity  being  minimized  being  the  N 
stage  local  cost  and  the  (N-1) -stage  optimal  cost  for  the 
resulting  state  or  stage  (N-1) . 


The  decision  at  stage  N  which  yields  the  minimum  total 
N-stage  cost  is  then  the  Nth  component  D (N)  of  the  N-stage 
optimal  policy  *D(N).  The  output  BOD  and  DO  resulting  from 
the  optimal  treatment  level  of  stage  N  allow  one  to  go 
directly  to  the  approprite  state  of  stage  (N-1)  and  select 
the  next  component  *D(N-l)  of  the  N-stage  optimal  policy. 
This  *  forward  pass'  continues  until  the  first  stage  is 
reached. 

At  this  point,  it  would  be  appropriate  to  mention  that 
dynamic  programming  is  simply  an  approach  to  be  used  in 
optimizing  serial  systems.  Dynamic  programming  is  unlike 
linear  programming  in  that  there  is  no  standard  mechanical 
procedure  for  structuring  and  solving  a  problem.  For 
example,  in  most  real  world  applications,  one  must  determine 
the  form  of  the  recursive  relationship.  Much  is  left  to  the 
individual  in  applying  the  technique  successfully. 
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CHAPTER  FOUR 


MATHEMATICAL  MODEL  OF  THE  SYSTEM 


4.1  Introduction 


In  this  chapter,  the  rationale,  concepts,  and 
mathematical  models  for  the  simulation  of  water  quality  in  a 
river  system  are  developed.  Emphasis  is  placed  upon  the 
theoretical  aspects. 

The  simulation  model  is  made  up  of  a  group  of 
mathematical  equations  which  are  linked  together  by  a 
progammed  logic.  The  purpose  of  the  model  is  to  generate  a 
reasonably  accurate  representation  of  the  stream  flow  and 
oxygen  balance  in  a  river  system.  The  model  determine  the 
mean  and  standard  deviation  from  the  historical  data  for 
each  waste  disposal  station.  The  historical  data  are  daily 
gage  records  from  gaging  stations  located  at  the  top  of  each 
reach.  The  model  simulates  quantity  of  waste  flow  into  the 
river  system  using  the  mean  and  standard  deviation  computed 
from  historical  data  for  each  reach1.  The  value  of  the 
simulation  model  depends  to  a  great  extent  upon  the  use  of 
high-speed  digital  computer  42  fast  and  accurate 


i  Described  in  detail  in  section  5.3.  ’Preparation  for 
Simulation  * . 
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computations  and  application  of  logic.  The  speed  of  the 
computer  not  only  provides  the  results  within  an  acceptable 
time  but  also  reduce  materially  the  costs  to  obtain  the 
results.  In  short,  the  computer  is  a  necessary  appurtenance 
in  the  use  of  the  simulation  model.  For  this  reason,  the 
development  of  models,  techniques  and  logic  is  made  for 
direct  application  of  computer  methods. 

A  system  in  which  there  are  several  waste  dischargers 
is  considered  to  be  divided  into  reaches,  with  continuously 
treated  or  untreated  waste  discharged  into  it.  The  reach  is 
defined  as  the  stretch  of  a  stream  between  two  dischargers. 
The  last  reach  encompasses  the  remainder  of  the  stream  from 
the  last  discharger  to  the  end  of  the  segment  under 
consideration.  The  optimal  allocation  of  artificial 
aeration  equipment  along  the  section  is  to  be  determined. 
This  is  done  by  starting  with  the  Streeper-Phelps  model 
(1925)  and  by  making  appropriate  simplifying  assumptions 
such  as  constant  river  coefficients,  steady-state  one 
dimensional  flow,  small  effect  of  second  stage  biochemical 
oxygen  demand,  negligible  dispersion  effect  and 
photosynthesis. 
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4 . 2  Theoretical  Model 

The  first  dissolved  oxygen  model  for  predicting  oxygen 
balance  in  a  flowing  stream  was  presented  by  Streeter  and 
Phelps.  Fcr  the  present  workr  the  simple  first-order  decay 
function  is  used  for  BOD  (4.  1)  and  the  Streeter-Phelps  model 
is  used  for  DC  {4.  2)  . 


dl 

dt 


(4.1) 


dc 

dt 


(Cs  -  Ct>  -  klLt 


(4.2) 


It  is  evident  that  any  particular  solution  to  equation 
4.1  is  independent  of  eguation  4.2,  whereas  the  solution  to 
equation  4.2  dependent  upon  the  solution  to  equation  4.1. 
That  is,  the  BCD  profile  is  independent  of  the  oxygen 
profile,  but  the  oxygen  profile  depends  upon  the  EOD 
profile . 

Combining  equations  4.1  and  4.2,  we  have 


dc 

dt 


k2  Cs  k2Ct  kl  L0  6 


“klt 


(4.3) 


With  artificial  aerator  at  point 
following  boundry  conditions: 


and  subjected  to  the 
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C  =  CQ  at  t=0 

~  ^i-l  ~  at  t  =  t^_^  (4.5) 

the  CO  concentration  at  any  time  t  is  given  by  the  following 
Streeter-Phelps  eguation. 


0  e  =  C  (e  -  1) - --(e  -1)  +  C  +V.e  2  1  V 

c  s  If  -L-  Oil 

2  kl 


(4.6) 


The  rate  of  oxygen  absorption,  U,  by  stream  water  depends  on 
numerous  interacting  factors,  such  as  water  temperature, 
flew  condition,  characterstics  of  water  and  type  of  air 
applied.  According  to  C’Ccnnor  and  Dobbins  (1956),  the  rate 
of  oxygen  supply  by  an  aerator  device  is  represented  by 


where  K  is  the  volumetric  absorption  rate  coefficient  of 

v 

oxygen.  This  depends  on  the  stream  temperature.  Morgan  and 

Eewtra  (1960)  carried  out  an  experiment  with  compressed  air 

to  evaluate  K  at  20°C  and  found  the  value  to  be  34.88  for 

v 

their  specific  eguipment.  Then  the  rate  of  artificial 
aeration  absorption  term  of  Morgan  and  Bewtra  can  be 
represented  as 

(t-20) 

U  =  34.88  (  1.024)  (C  -  C  ) 

s  t 

The  following  equation  is  obtained  by  just  rearranging  the 
eguation  4.6. 
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k2t 


(Cs-Ct)e  '  (Cs'C0)+ 


klL0  ,  (k2  kl)t: 


k2"kl 


(e 


-  1)  -  U.  V.e 
x  1 


k2ti-l 


(4.7) 


By  introducing  tie  DO  deficit  instead  of  DO 
concentration,  eguation  (4.7)  becomes: 


'k2t  klL0  "V  "k?t)  -k9(t-t.  ,) 

D.=  D„e  +  ------ (e  1  -  e  2  )  -  U  V  e  2  1-1 

k  -k  1  i 

2  kl 


0 


(4.8) 


Ihis  is  the  governing  equation  for  DO  sag  profiie  along 
a  stream  with  a  point  aeration  source  at  (i-l) 


for  dynamic  programming  formulation  convenience,  the 

river  segment  under  consideration  is  divided  into  N  reaches, 
ttl 

the  i  reach  stretches  from  ti_1  to  t ±.  The  length  of 

reach  equals  (t^  -  t  i_1)  =  (days)  .  Assume  that  aerators 

A  q,  A-^ ,  A^,  ■  ......  ,  ar  e  1  o c at ed  at  0,  1,  2,...., 

i, . . . .  N- I ,  the  top  of  each  of  N  reaches.  For  simplicity, 

cacacities  of  these  aerators  are  assumed  to  he  equal  to  the 

increment  in  level  of  DC  needed  at  chose  points,  they  are 

V  ,  V  . ,V  ....V  .  The  Figure  4.1  gives  the  schematic 

12  i  N 

representation  of  the  system. 


' 
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4.  1 

Schematic  representation  of 

the  system. 

The  DO 

deficit  at  any  point  ? 

within  the 

reach  i  can  be 

computed  by  equation  (4.9)  where  Di_1  is  the  initial  DO 
deficit  for  this  reach, 
let  p  >  0  and  p  <  Xi  days. 


„  -  n  / 2P  klLi-l  ,  'klP  _k2P.  „  _  "V 

Di-l+p-  Di-le  +  (e  *  e  >  -  Ui  Vi  e 

k2-  kx 


(4.9) 


Where 


<  p  <  X. 


At  the  end  point  of  this  reach/  P  =  X  ,  the  DO  deficit 


becomes  Di  . 


-k_X.  knL.  ,  -k  X.  -k0X. 

r,  /  T\  TT  TT  N  2  1  ,  1  1-1  ,  1  1  2  lv 

D.  =  (D.  -  V.  U.)  e  + - (e  -  e  ) 

1  1-1  1  1  k  -  k 

k2 


(4.10) 


• 
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Again  from  reach  to  reach  the  BOD  changes.  The  BCD 
concentration  at  the  top  of  each  reach  depends  on  the 
guantity  of  waste  added  at  that  point.  Accordingly,  the  DO 
deficit  varies.  The  BOD  concentration  at  the  end  of  the 
reach  is  given  by 


~k. X.  0  -k1X. 

Li-l=Li-2e  =Li-3e  e 


-(i-l)X.k 

=  L0  e  (4.11) 

Where  1 n  is  the  initial  BOD  (BOD  )  level  plus  the  BOD 
u  in 

(BODdis  )  at  the  top  of  the  reach.  Substituting 

equation  (4.11)  into  equation  (4.10)  gives 


-(i-l)kX, 

-k  X  k  L  e  1  -k  X. 

D,  =  (D._1-  U.V.)  e  2  1  +  -i-2 -  (e  1  1 


-  e 


■k2Xl> 


k2  -  k 


(4.12) 


This  equation  gives  the  DO  deficit  at  the  end  of  the 
th 

i  reach.  The  deficit  should  not  be  less  than  the  standard 
value  specified  by  the  authorities.  In  order  to  maintain 
the  value  of  DO  greater  than  or  equal  to  the  specified  value 
(5  mg/1) ,  artificial  aeration  should  be  provided  at  the  top 
of  this  reach.  We  must  provide  aeration  to  the  stream  in 
such  a  way  that  the  specified  law  value  (C  >  5mg/l)  should 
be  maintained  for  all  parts  of  the  stream  in  that  reach.  No 
distinction  is  made  between  day  time  and  night  time  operation. 


hba 
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The  capacity  of  aerator  located  at  A  can  be  computed 
from  eguation  (4.12),  provided  that  the  value  of  D  and  D 
are  specified,  and  all  other  constant  and  coefficients  are 
given,  that 


kiLoe 


-(i-l)k.X. 

1  i 


(k  -k  )X 

(e  1  1)  +  D._1- 


D.e 

l 


k0X. 
2  i 


(4.13) 


Ihe  above  equation  gives  the  amount  of  aeration  needed 
to  keep  the  DO  level  above  the  specified  value  in  that 
reach.  Otherwards  the  capacity  of  aerators  at  the  top  of 
the  i  reach. 


5 


■  Aim 
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4.3  The  Objective  Function  and  Constraints 

Ihe  second  objective  is  to  define  the  objective 
function  to  achieve  the  artificial  aeration  at  minimum  cost. 
For  this  aeration  problem  the  mathematical  form  of  the 
objective  function  is  written  as 

N 

z  =  l  wi  (ci  ■  +  w2  v±  (4-i4) 

Where  Wx  and  «2  are  weighting  factors  («1  >  0,  W  >0  and 
W  ±  +  W  2  =  1)  and  N  is  the  number  of  aerators  used  in  the 
system . 

The  first  term  represents  the  weighted  sum  of  the 
squares  of  the  specified  law  value  of  DO,  and  the  actual  DO 
level  C  The  objective  function  should  taken  on  a  zero 

value  whenever  the  specified  law  value  is  reached  for  all 
time.  The  corresponding  value  cf  would  also  be  zero. 
Whenever  (C  -  )  =0,  the  objective  function  should  always 

be  positive.  If  Ci  <  C  damaging  cost  for  the  environment 
is  incurred  and  if  C.  >  C0  ,  unnecessarily  providing  aeration 
cost  is  imposed.  The  costs  can  be  either  the  time  cost  in 
dollars  per  unit  DO  level,  if  any  estimation  method  is 
available  in  the  literature  or  the  nominal  costs  used  in  a 
typical  objective  function  in  the  control  problem. 
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The  use  of  weighting  factors  in  the  objective  eguation 

provides  the  planner  a  planning  flexibility.  For  instance, 

if  the  designer’s  main  concern  is  to  satisfy  the  stream 

standard,  then  he  can  select  in  such  a  way  the  ratio 

a  ^ar9e  numbe£  ob  if  the  main  concern  is  the  cost 

for  providing  artificial  aeration  and  the  water  quality  is 

permitted  to  deteriorate  at  certain  points  along  the  stream, 

then  he  can  choose  W  such  that  W  /W  is  large.  Hence  by 

2  2  1 

varying  different  combinations  of  the  values  of  W  and  W2 
provides  the  designer  a  designing  flexibility. 

since  (C.-  Cz)  =  (Cs  -  C£)  -  (Cg  -  C.  )  =  D  %  -  D.  ,  eguation 
(4.  14)  can  be  expressed  in  terms  of  DO  deficit,  that  is 


Minimize  Z(D,V) 
v 


N 

I 

i=l 


wi  (V  Di)2  + 


(4.15) 


Ey  employing  eguation  (4.  15)  the  mathematical  model  of  our 
problem  can  he  formulated  as  follows: 


Minimize  Z(D,V) 
v 


wi  <V  V2  +  v2 

i=l 


(4.16) 


Subjected  to: 


I.  The  limitation  on  total  available  aeration  capacity 


9  *  ' 


aaolioi  as  bssfcXiii.ici  »d  a  so  aaXdciq 
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N 

Z  V.  <  Y 

•  i  l  — 


(4.17) 


where  Y  is  the  available  aeration  capacity. 

2.  Physical  limitation  on  DO  deficit  and  aeration  capacity. 

0  <  v±  <  Cs  (i  =  2,  . .  N)  (4.  18) 

o  <  v.  <  cg  (i  =  lr  2,  . ,  N)  (4.  19) 

Ey  introducing  a  Lagrange  multiplier  €  (€  >  0)  ,  which 

has  the  significance  of  a  price  (Bellman,  1962) ,  the 

modified  objective  function  eguation  (4.20)  with  constraints 

(4.2  1)  and  (4.22)  below  is  equivalent  to  the  original 

problem  (eguation  4.16,  4.17,  4.18  and  4.19)  (Everett, 

1963)  ,  where  Y  eguales  z  V  found  in  the  modified  problem. 

i 

N  2  2^ 

Minimize  Z(D,V)  =  Z  W_ (D0  -  D.)  +  W~V7  +  €  E  V.  (4.20) 

v  i=l  1  £  1  2  1  i=1  i 

Subjected  to: 

0  <  D  .  <  C  (h.2  I) 

i  s 

0  <  V  .  <  C  (4.22) 

l  s 

Since  the  optimal  policy  of  the  modified  problem  is  a 
function  of  €  ,  We  can  find  a  proper  value  of  e  by 
interpolation  to  force  the  optimal  policy  to  satisfy  the 
level  of  Y  desired  in  constraint  (4. 17). 
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4.4  Derivation  of  the  Securrence  Relation 

Once  the  objective  function  has  been  formulated,  the 
next  step  in  the  mathematical  formulation  centers  around  in 
developing  the  recurrence  relation. 

As  defined  before,  the  reach  of  interest  is  divided 
into  N  stages  where  stage  i  corresponds  to  the  end  of  reach 
i.  Let  D#  be  the  state  variable  at  stage  i,  it  is  the  DO 


1 


deficit  at  point  t  and  V  be  the  decision  variable  at  stage 

i  i 

i,  it  is  the  DC  increment  needed  at  point  i  -  1,  when  D 
and  E  ^  are  specified. 

The  stage  transformation  eguation  T  can  be  derived  from 
equation  (4.  10) 


or 


-(i-  l)Xik1 


(k0-k  )X 


Di_r  D.e 


+  V.- 

X 


2  r  i 


-  i) 


(4.24) 


The  recursive  formula  is: 


With  the  following  conditions 


0  <  VN  <  c 


S 


(4.26) 


9  io  f»d: 
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0  <  E  <  C 
N  s 


(4.27) 


The  objective  function  for  the  first  stage  becomes: 


Z-j  (D-j  ,6)  =  Min  W1(D1  -  D^2  +  V?  +  6  Vi 


v. 

i 


where 


0  <  V.  <  C 
i  s 

0  <  D.  <  C 
1  S 


(4.28) 


(4.29) 

(4.30) 


In  equation  (4.25)  the  subscript  N  represents  the 
number  of  stages  remaining  in  the  process;  the  vector 
describes  the  state  of  the  process  with  N  stages  remaining; 
the  function  Z^(D^,g  )  is  the  yield  from  an  N  stage  process 
starting  at  state  S  and  employing  an  optimal  policy;  is 
the  class  of  all  admissible  policies  or  decisions; 

Z  i (D  i ,  6 )  is  the  yield  during  the  first  stage  of  the  N 
stage  process  and  depends  upon  decision  .  Mathematically, 
we  have  said,  if  one  knows  the  immediate  effect  of  a 
decision  over  the  next  time  interval  and  the  long-range 
effect  from  the  end  of  that  time  interval  until  the  end  of 
the  process,  then  one  knows  the  effect  of  that  decision  from 
the  present  until  the  end  of  the  process.  This  is  precisely 
the  total  return  over  the  process. 


' 


CHAPTER  FIVE 


APPLICATION  TO  NORTH  SASKATCHEWAN  RIVER 


5.  1  Introduction 


The  simulation  model,  including  the  data  generator,  was 
applied  to  the  North  Saskatchewan  River  near  Edmonton, 
Alberta,  Canada.  This  river  holds  an  unique  position  among 
the  major  rivers  on  the  North  American  continent  in  that,  as 
a  major  source  of  water  supply  and  as  a  receiver  of 
municipal  and  industrial  waste  water  (Figure  5. I) ,  it  is 
completely  frozen  over  during  the  winter  for  a  periods 
averaging  five  months.  During  this  period  of  ice  cover, 
minimum  discharges  also  occur  in  the  river. 

During  the  winter  of  1955-56,  the  DO  level  100  miles 
downstream  from  Edmonton  was  zero  (Figure  5.2).  This  led  to 
an  attempt  to  artificially  aerate  the  river  (Department  of 
Public  health,  1956).  With  the  cooperation  of  the  city  of 
Edmonton  and  the  Provincial  Department  of  Public  Works,  the 
Sanitary  Engineering  Division  of  the  Department  of  Public 
Health,  set  up  an  air  compressor  system  near  Fedwater, 
Alberta  in  an  attempt  to  artificially  raise  the  dissolved 
oxygen  levels  in  the  river.  The  operation  continued  for 
eighteen  days,  at  which  time,  melting  conditions  made  it 
necessary  to  remove  equipment  from  the  ice.  During  the 
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FIGURE  5.1:  Map  of  North  Saskatchewan  River  in  Alberta 
and  Waste  Outfalls  on  the  River  Banks. 
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continuous  operation,  a  raise  of  only  0.25  mg/1  in  dissolved 
oxygen  (DO)  level  was  observed  immediately  downstream  from 
the  aeration  system.  No  increase  in  DO  was  observed  further 
downstream.  This  led  to  the  introduction,  by  the  Alberta 
Division  of  Environmental  Health  Services,  of  stringent 
control  of  all  wastes  discharged  to  the  river.  As  a  result 
of  this  pollution  abatement  program,  improved  dissolved 
oxygen  levels  have  been  maintained  in  the  river.  This  has 
been  brought  about  by  the  construction  of  the  primary 
treatment  section  of  a  new  sewage  treatment  plant,  augmented 
in  1957  by  the  use  of  secondary  treatment  facilities  at  this 
plant.  In  1961,  Brazeau  Dam  was  constructed  to  provide 
storage  for  winter  release  to  improve  water  quality  during 
winter  months.  In  1965,  sewage  lagoons  were  constructed  for 
winter  storage  for  meat  packing  plant  wastes  and  some 
domestic  sewage.  The  use  of  these  lagoons,  sewage  treatment 
plant,  and  low  flow  augmentation  dam  have  further  reduced 
the  waste  loading  on  the  river. 


.  i 
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5 • 2  River  Flow  an  d  Its  Usage 

The  North  Saskatchewan  River  is  characterized  not  only 
by  its  relatively  long  period  of  ice  cover  but  also  by  the 
wide  variation  between  maximum  and  minimum  daily  flows. 

This  is  shown  in  Figure  5.3a,  5.3b  and  Table  5.1.  Due  to 
this  broad  variation  in  flow  and  continuous  constant  daily 
waste  discharge  into  the  river  system  makes  it  difficult  to 
maintain  the  guality  of  water  throughout  the  year.  This 
uncertainity  in  river  condition  is  an  ideal  problem  for 
investigation  by  simulation. 

Further,  this  river  serves  as  a  source  of  water  supply 
for  domestic  consumption,  agricultural  and  manufacturing 
processes  and  as  a  receiver  and  carrier  of  domestic  and 
industrial  wastes.  Table  5.2  lists  the  major  users  of  the 
river  in  the  province  of  Alberta. 

Because  of  its  large  population  and  industrial  complex, 
the  city  of  Edmonton  is  the  source  of  major  pollution  load 
to  the  North  Saskatchewan  River  in  Alberta.  The  main  source 
is  the  domestic  sewage  effluents  from  the  main  and  No. 3 
sewage  treatment  plants. 

A  series  of  sewage  lagoons  are  located  16  miles 
downsteram  from  the  low  level  bridge  at  Edmonton.  These 
lagoons  receive  domestic  and  meat  packing  house  wastes  from 
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jURE  5.3: 


AVERAGE  MONTHLY  FLOW 
NORTH  SASKATCHEWAN  RIVER 


DATA  FROM  FEDERAL  DEPARTMENT  OF  WATER  RESOURCES 


Maximum  mean  daily  *  164,  OCO  c.f.s.  (June  28,  1915) 
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the  northeast  section  of  the  city  and  domestic  wastes  from 
the  hamlet  of  Sherwood  Park,  is  one  of  the  major  users  of 
this  river. 

The  Edmonton  industrial  complex,  consisting  mainly  of 
the  petroleum  and  chemical  industries,  use  the  North 
Saskatchewan  River  as  a  source  for  water  supply,  cooling 
process  and  as  a  receiver  of  industrial  waste  effluents.  In 
the  past,  these  industrial  wastes  have  been  the  main  cause 
of  taste  and  odor  problems  found  downstream  from  Edmonton. 

Recreational  use  of  the  North  Saskatchewan  River  in 
Alberta  is  restricted  to  boating  and  some  sport  fishing. 
There  is  no  commercial  fishing  undertaken  in  this  river  in 


Alberta. 


. 
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TABLE  5. 1 


NORTH  SASKATCHEWAN  RIVER  DISCHARGES 


(Information  from  the  Federal  Water  Resources  Branch) 


YEAR 

Maximum  in¬ 
stantaneous 
discharge 

c .  f .  s. 

191  1 

1912 

1913 

1914 

19  15 

204500 

1916 

1917 

1918 

19  19 

1920 

61600 

1921 

2  74  00 

1922 

28600 

1923 

99600 

1924 

27600 

1925 

77000 

1926 

1927 

1928 

1929 

1930 

23900 

193  1 

1932 

1933 

1934 

1935 

1936 

1937 

1938 

1939 

1940 

Maximum  Date  of 
daily  maximum 

discharge  daily 
c .  f .  s.  discharge 


5  14  42 

Jul. 

3 

74100 

J  ul . 

10 

32600 

Aug. 

15 

6  1740 

J  un. 

9 

164000 

Jun. 

29 

58800 

Jun. 

22 

65597 

May 

18 

35347 

Jun. 

16 

19885 

Jun. 

24 

57220 

May 

10 

24888 

May 

23 

25760 

Aug. 

18 

84100 

Jun. 

25 

27500 

Jul. 

5-6 

75800 

Aug. 

18 

58700 

Sep. 

4 

40400 

Jun. 

29 

6  12  00 

Jul. 

7 

38100 

Jun. 

5 

23700 

Jul. 

17 

39200 

Jul. 

2 

66000 

Jun. 

4 

34400 

Jun. 

19 

28  100 

Jun . 

1 

46300 

Jul. 

1  1 

40400 

Apr. 

19 

3  1500 

Jul. 

1  7 

40000 

J  ul. 

4 

30200 

Jun. 

28 

35700 

Apr. 

18 

Minimum 

Date  of 

daily 

minimum 

discharge 

daily 

c.  f .  s. 

discharge 

984 

Nov.  12 

1062 

Mar.  25 

12  10 

AP.  6 

650 

Dec. 24-26 

700 

Dec.  15 

950 

Mar.  4 

1  100 

Feb.  22 

960 

D. 1 2-Fe. 2  0 

688 

Mar.  4 

895 

Dec.  4 

800 

Dec.  23 

380 

Nov.  25 

540 

D.  12, J.  1  1 

760 

Jan.  24 

538 

Feb.  26 

1  140 

Jan.  19-31 

1290 

Dec .  16-18 

1330 

Nov.  13 

965 

J. I9-F.2 

952 

J.  1 6-F.  14 

700 

Mar.  16 

865 

Mar.  14 

796 

Dec.  13 

900 

Mar.  2 

406 

Jan .  3 

485 

Apr .  4 

480 

D. 23, J. 3 

682 

Feb.  4,14 

1030 

Mar.  3-5 

220 

Jan .  1 
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YEAR 

194  1 

1942 

1943 

1944 

1945 

1946 

1947 

1948 

1949 

1950 

1951 

1952 

1953 

1954 

1955 

1956 

1957 

1958 

1959 

1960 

196  1 

1962 

1963 

1964 

1965 

1966 

1967 

1968 

1969 

1970 
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(Cont.  TABLE  5.  1) 


Maximum  in- 

Maximum 

Date  of 

Minimum 

Date  of 

stantaneous 

daily 

maximum 

daily 

minimum 

discharge 

discharge 

daily 

discharge 

daily 

c .  f.  s . 

c .  f .  s . 

discharge 

c.  f .  s. 

discharge 

26720 

Jun.  28 

548 

Nov.  17 

43960 

42250 

Jul.  14 

583 

Nov.  26 

44020 

Apr.  12 

1080 

Feb.  10-13 

125900 

12  1970 

Jun.  16 

551 

Nov.  22 

25060 

24300 

Jun .  1 

724 

Nov.  30 

44760 

Jun.  24 

1300 

Feb.  8 

28600 

Jun.  13 

602 

Nov.  27 

66620 

6  54  40 

May  25 

1  140 

Feb.  12 

32680 

Jul.  22 

730 

Dec .  1 4 

53720 

50330 

Jun.  17 

430 

Dec .  1 3 

40980 

39020 

May  3 

624 

Nov.  25 

132000 

125000 

Jun.  25 

1030 

Dec.  23 

45800 

44900 

Jun.  5 

652 

Nov.  29 

118400 

106600 

Jun.  8 

833 

Apr .  8 

32020 

30380 

Jun.  15 

1040 

Jan .  5 

26600 

25460 

Jun.  7 

580 

Nov.  20 

23380 

21780 

Jun.  11 

506 

Dec.  13 

52  130 

49890 

Jul.  1 

1  180 

Jan.  7-8 

51740 

46  140 

Jan.  29 

598 

Dec .  1 4 

388  10 

36830 

Jul.  3 

640 

Nov.  20 

30  100 

272  10 

Jul.  31 

700 

Nov.  30 

28500 

27000 

Aug.  6 

575 

Nov.  28 

39900 

37  100 

Jul.  18 

1330 

Jan .  2 

49700 

4  7600 

Jun.  2  1 

1350 

Nov .  2  5 

95300 

9  1600 

Jun.  29 

1070 

Dec.  1 

57800 

Jul.  6 

2280 

Dec .  1  5 

5  1200 

Jul.  9 

1500 

Dec  •  7 

50769 

43580 

Aug . 

1820 

Dec.  20 

57209 

Jun.  18 

1097 

Nov.  27 

93124 

80563 

Jul.  10 

2100 

Nov .  8 

• 
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TABLE  5.2 

M A JOE  USERS  OF  THE  NORTH  SASKATCHEWAN  RIVER 

(Information  supplied  by  the  Alberta 
Department  of  Public  Health) 


Users 

Purpose  of  Use 

Town  of  Rocky  mountain 
house 

Water  supply  &  sewage  disposal 

Town  of  Drayton  Valley 

Water  supply  &  sewage  disposal 

Town  of  Devon 

Water  supply  &  sewage  disposal 

Imperial  oil  -  Devon 

Industrial  waste  disposal. 

City  of  Edmonton 

Water  supply,  sewage  disposal, 
snow  dumping  &  thermal 
discharge  from  power  plant. 

Edmonton  Industries: 


Canadian  Industries  Ltd. 
TEXACO 

Building  Products  Ltd. 

S  &  L  oil  Ltd., 

Union  Carbide, 

Gulf  Oil  Canada  Ltd., 
Uniroyal  Ltd. , 

Imperial  Oil  Ltd. 

McColl  Frontenac  Oil  Co. 
British  American  Oil  Co. 
Chemcell  Ltd. 

Water  supply  and 

Industrial  waste  disposal. 

County  of  Strathcona 

Municipal  Sewage  disposal. 

City  of  Edmonton  Lagoons: 


Alberta  Hospital  -  Oliver 
Sherwood  park 

Edmonton  meat  packing 
plants. 

Domestic  &  Industrial 

Waste  storage  for  disposal 
in  summer  months. 

Town  of  Fort  Saskatchewan 

Municipal  sewage  disposal 

Sherritt  Gordon  Mines  Ltd. 
Dow  Chemical  Ltd. 

Water  supply  and 
waste  water  disposal. 

_  *:  ,  ~ 
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(cont.  TABLE  5.2) 


Users 

Town  of  Redwater 

Imperial  Oil  -  Redwater 

Western  Chemicals  Ltd. 
Duverney . 

Town  of  Elk  Point 

Canadian  salt  Co.  Ltd. 
Lindberg. 


Purpose  of  Use 

Water  supply  and 
waste  water  disposal. 

Industrial  waste  disposal. 

Water  supply  and 
Waste  water  disposal. 

water  supply  and 
waste  water  disposal. 

Water  supply  and 
waste  water  disposal. 


. 

. 

. 
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5.3  Preparation  for  Si  mutation 

The  river  flow  simulation  requires  a  long  and 
continuous  record  of  several  years  of  river  flow  data. 
Further,  the  river  system  under  consideration  may  be  stated 
as  follows  for  dynamic  programming  convenience. 

The  first  step  is  to  discretize  the  stream  system  by 
dividing  it  into  reaches.  The  system  is  consider  to  be  made 
up  of  a  number  of  connected  segments  in  which  each  segment 
has  constant  properties.  Peach  nodes  are  therefore 
established  at  every  major  waste  discharge  point  or 
significant  change  in  the  system.  Significant  changes  are 
the  tributary  flows,  changes  in  flow  conditions,  etc. 

The  segment  of  the  river  considered  in  this  work  is 

about  200  miles  long.  This  segment  is  divided  into  24  (N ) 

th 

reaches  as  shown  in  Figure  5.4.  The  i  reach  stretches  from 
t  i-ito  tf.  The  time  of  travel  in  each  reach  depends  on  the 
velocity  and  other  hydraulic  constants.  The  time  or  the 
length  of  a  recah  is  equal  to  (%_i-  t ±  days,  time 

expressed  in  length.  Assume  that  aerators 
A  o  /  ^2  ,•••••,  Aj^_^  are  located  at  0,I,2,....,N<~I,  t  he 

top  of  each  of  N  reaches.  For  simplicity  the  capacity  of 
these  aerators  are  assumed  to  be  equal  to  the  increment  in 
the  level  of  DO  needed  at  those  points,  they  are  V1 , 

v2  . . .  %  • 
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DISTANCE  DISPOSAL  LENGTH  OF 

IN  MILES  STATIONS 

MILES 

0.0  EDMONTON 

REACHES  TYPE  OF  WASTE  DISTANCE 

DISCHARGED  IN  DAYS 

DAYS 

DOMES.  +  INDUST.  WASTE  0  tiayq 

7.6  miles 

1.6  CLOVER  BAR 

'  .63  DAYS 

DOMES.  +  INDUST.  WASTE  .63 

8.0 

13.6  LAGOONS 

:  .67 

•  • 

INDUSTRIAL  WASTE  1.3 

7.3 

22.9  FORT  SASKATCHWAN 

:  .61 

:  DOMES.  +  INDUST.  WASTE  1.91 

28.1  SHERRITT  GORDON  5,2 

•  43 

*  DOMES.  +  INDUST.  WASTE  2.U 

32.5  STURGEON  RIVER  4,4 

•  30 

:  RIVER  FLOW  2.64 

6.2 

38. 7  VINCA 

:  .si 

:  DOMES.  +  INDUST.  WASTE  3. IS 

41.9  TCWN  OF  REDWATER  TT2 

;  .27  DIMES.  +  INDUST.  WASTE  3. A? 

4664  ELDORENA 

:  .31  DOMESTIC  WASTE  3.77 

9.6 

56.0  WASKATENAU 

• 

:  .so 

:  DOMES.  +  INDUST.  WASTE  4. S3 

7.3 

63.3  WARSPITE 

:  .6i 

:  DOMESTIC  WASTE  5.14 

9.9 

73.2  PAKAN 

:  .83 

DOMESTIC  WASTE  5.97 

8.0 

81.2  WHITE-EARTH  CREEK 

.63 

TRIBUTARY  6.6 

8.3 

89.5  SHANDRO 

.70 

DOMESTIC  WASTE  7.3 

9.6 

99.1  DESJARLAIS 

.80 

DOMESTIC  WASTE  8.1 

17.0 

116. 1  DUVERNAY 

1.42 

DOMES.  +  INDUST.  WASTE  9.52 

10.8 

126.9  BEAUVALLON 

.83 

DOMESTIC  WASTE  10.35 

8*  0 

134.9  MYRNAM  * 

.63 

DOMESTIC  WASTE  11.02 

9*6 

144.5  HOPKINS 

.80 

DOMESTIC  WASTE  11.82 

6.1  : 

150.6  ELK  POINT  I 

DOMES.  +  INDUST.  WASTE  12.32 

10.6  : 

161.2  LINDBERGH  ! 

.88 

DOMES.  +  INDUST.  WASTE  13.2 

8.9  | 

170.1  HEINSBEROH  : 

.75 

DOMESTIC  WASTE  13.95 

U,  | 

181.6  LEA  PARK  1 

.95 

DOMESTIC  WASTE  14.9 

ua  j 

194.7  FORBESVILLE 

1.10 

DOMESTIC  WASTE  16.0 

6.3 

201.0  LLOYDMINSTER  : 

J’  ALBERTA-SASK.  BORDER  16.5 

FIGURE  5.4:  Schematic  Illustrating  the  North  Saskatchewan  River 
Division  of  Reaches,  Length  of  Reaches,  Type  of  Waste 
Discharged,  Distance  Downstream  from  Edmonton. 
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There  is  no  real  reason  to  limit  the  number  of  reaches 
except  to  meet  a  constraint  on  machine  storage.  The  program 
coding  in  this  simulation  model  is,  for  this  reason,  limited 
to  24  reaches  which  should  be  adequate  for  most  quality 
control  studies. 

The  next  step  is  the  numbering  of  the  reaches. 

Although  an  initial  computational  sequence  is  included  in 
the  model,  and  numbering  may  be  in  any  pattern.  For  dynamic 
programming  application  it  is  most  convenient  to  start 
numbering  at  the  downstream  end  and  number  consecuti vely. 
Then,  continue  numbering  upstream  on  each  reach,  starting 
with  the  downstream  most  reach,  until  all  reaches  are 
numbered.  The  numbers  should  be  consecutive,  I,  2, 

3,  -  N  for  N  reaches  in  the  system. 

The  basis  for  this  simulation  of  the  stream  flows  is  an 
historical  record  of  fifty  years  daily  stream  flow  data.  In 
the  simulation,  river  flow  data  are  generated  in  a  manner 
such  that  the  statistical  properties  of  the  historical  data 
are  preserved.  It  is  desirable  to  have  a  long  and 
continuous  record  of  several  years  of  data.  The  primary 
source  of  streamflow  data  is  from  the  Federal  Water 
Resources  Department,  Edmonton  branch.  This  data  is  stored 
on  magnetic  tape  and  is  processed  further  for  use  in 
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The  model  used  in  this  work  was  developed  specifically 
for  generating  flows  at  each  stations.  The  historical  data 
are  daily  gage  records  from  gaging  stations  in  the  river 
basin,  each  covering  an  identical  period.  The  data  are 
checked,  missing  data  are  filled  and  monthly  averages  are 
computed  by  an  auxiliary  program  called  CHKDATA 
(Appendix-A)  . 

The  transformed  historical  monthly  gage  data  are 
analyzed  statistically  to  determine  the  mean  and  standard 
deviation.  This  mean  and  standard  deviation  computed  from 
historical  data  for  waste  disposal  station  i  is  used  to 
simulate  the  rate  of  waste  flow  into  the  system  at  stage  i 

(i  =  I,  2,  - ,  N=24).  The  river  flows  and  waste  flows 

from  each  stations  were  generated  by  making  a  linear 
transformation.  The  objective  is  to  maintain  the 
statistical  appearance  of  the  historical  data  in  generating 
the  synthetic  data. 
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5.4  Simulation  Model  Operation 

The  simulation  model  is  actually  two  simulation  models. 
The  first  model,  CHKDATA,  is  designed  to  read  daily  raw 
streamflow  data.  It  then  computes  monthly  averages  for  each 
reach.  The  monthly  averages  computed  by  CHKDATA  will  be  an 
input  to  the  second  model,  RIPSIM.  This  model  simulates  the 
riverflow  and  quantity  of  waste  discharged  at  each  reach. 
Then  it  calls  the  subroutine  OPTMOD,  to  optimize  the  amount 
of  aeration  required  at  the  top  of  ecah  reach  to  maintain 
water  quality.  The  two  programs  can  be  run  operated 
together  requires  considerable  computer  capacity  (42500 
bytes) .  Generally,  it  is  convenient  to  run  the  CHKDATA 
separately,  storing  the  output  on  magnetic  tape  (or  disk). 
This  output  then  becomes  input  to  the  simulation  program. 

The  simulation  program  requires  inputs  of  gage  data; 
reach  designations,  location,  length  of  each  reach,  initial 
BOD  and  DO  deficit,  minimum  DO  value  and  the  rate  of  waste 
flow  into  each  reach.  The  program  establishes  the 
computational  sequence  and  sets  up  the  weight  matrix  by 
computing  the  aeration  cost  for  the  discretized  DO  values. 
Simulation  begins  at  the  downstream  end  and  proceeds,  reach 
by  reach,  upstream,  computing  the  BOD  concentrations,  DO 
deficit,  at  the  upstream  and  downstream  ends  of  each  reach 
and  checks  for  a  deficit  stationary  point  within  the  reach 
under  consideration.  The  maximum  value  of  the  deficit  in 
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that  reach  is  subtracted  from  the  DO  saturation  value  (Cs) 
to  obtain  the  minimum  DO  concentration,  which  is  compared  to 
the  standard.  A  value  less  than  the  standard  results  in  an 
indicated  violation.  The  program  outputs  the  violation  DO 
value,  reach  number  and  the  time,  the  violation  occured.  To 
avoid  this  violation  the  optimization  program,  CPTMOD, 
compute  the  amount  of  aeration  required  at  the  top  of  that 
reach. 


In  the  actual  program  constraint  (4.  IS)  was  checked  at 
each  stage  by  two  IF  statements  to  force  the  aeration 
capacity  V(i,j,k)5  to  satisfy  the  constraint  (4.19)  in  the 
following  ways: 

If  V(i,j,k)  was  larger  than  Cs,  then  set  it  equal  to 
Cs,  if  it  is  less  than  zero,  then  a  very  large  cost  is 
assigned.  This  negative  aeration  capacity  will  definitely 
not  constitute  an  optimal  policy  due  to  its  very  large  cost. 
Constraint  (4.18)  can  be  relaxed  by  defining  a  feasible 
region  for  the  state  variable  (0  to  9mg/l) .  In  this 
feasible  region,  the  state  variable  will  be  discretized  into 
M  (M=  1  1)  values  with  Cs  as  its  maximum  and  zero  as  its 
minimum  value. 


5  Defined  in  Appendix-A,  section  A2.4  dictionary  of 
Variables1  . 


. 


. 


. 

. 


.  • 


.  • 


73 


To  start  the  algorithm,  the  number  of  aerators  are 
assumed  to  be  equal  to  the  number  of  stages  N.  At  each 
stage  a  decision  regarding  the  amount  of  aeration  must  be 
made.  The  amount  of  aeration  vary  from  0  to  9mg/l  in  this 
case. 


At  any  stage  i  each  one  of  the  M  discretized  states 
D(i,k)  is  coming  from  any  one  of  the  M  discretized  states 
D(i"Uj)  of  stage  i- 1  by  a  corresponding  aeration  with 
capacity  V(i,j,k)  to  each  K  and  subsequently  there  are  M 
total  cost  values  Z(i,j,k)  for  each  K.  Among  the  M  cost 
values  for  a  particular  K*,  the  minimum  cost  and  the  state 
position  j*  which  gives  the  minimum  total  cost  for  this 
particular  K*  are  stored.  We  store  such  M  minimum  cost 

values  corresponding  to  K*  =  1,  2,  3,  - ,  M  at  stage  i. 

At  the  final  stage  N,  we  pick  up  the  minimum  of  those  M 
minimum  cost  values  from  a  particular  position  *K  (K*  and  *K 
are  different,  where  *K  is  the  minimum  of  K*  values)  which 
is  the  minimum  total  cost  of  the  whole  system. 

By  tracing  forward  the  corresponding  state  positions 
which  eventually  give  the  final  state  position  *K  at  every 
next  stage.  This  gives  the  optimal  policy  for  that  selected 
e  value.  If  this  optimal  policy  satisfy  the  constraint 
(4.17),  then  it  is  the  optimal  policy  for  the  original 
problem,  equations  (4.  12)  and  (4.  16)  to  (4.  19)  .  If  not,  the 
value  of  6  is  modified  and  the  whole  procedure  is  repeated. 
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A  new  unconstrained  minimization  system,  equations  (4.20)  to 
(4.22)  is  formed  whenever  a  new  nonnegative  €  is  used.  If 
an  unconstrained  minimum  of  this  new  system  can  be  found  to 
be  Z*  and  the  total  aeration  capacity  used  in  achieving  this 
optimal  solution  to  be  Y*,  then  by  the  Everettfs  theorem 
(Everett  III  1963)  this  solution  is  a  solution  to  the 
original  problem,  equations  (4.  12)  and  (4.  16)  to  (4.  19)  , 
except  Y*  is  used  instead  of  Y  in  equation  (4.17) .  Hence  by 
sweeping  through  a  set  of  €  values  to  solve  a  spectrum  of 
problems,  total  available  aeration  capacities  and  their 
corresponding  optimal  policies  are  produced  in  the  course  of 
solution . 

Instead  of  assigning  a  fixed  value  to  the  total 
available  aeration  capacity  Y,  17  different  values  for  € 
were  assigned  over  the  range  between  zero  and  one  hundred 
and  then  solved  each  of  the  17  unconstrained  problems. 

The  entire  process  can  be  considered  as  a  learning 
device.  Each  value  of  6  yields  a  certain  amount  of 
information  about  the  system  response.  This  information  is 
considered  for  the  feedback  to  the  simulator  and  the 
simulator  generates  further  action  and  the  cycle  is 
repeated.  This  process  is  continued  until  an  acceptable 
solution  is  obtained.  The  existance  of  optimal  value  for  an 
objective  function  is  assured  by  the  algorithm,  since  all 
possible  ways  have  been  evaluated  and  compared  by  the 


. 
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algorithm.  For  this  application,  there  are  24  stages  with 
11  discretized  states  at  each  stage  and  the  dynamic 
programming  algorithm  find  the  optimum  value  out  of  all  II24 
possible  aeration  policies  in  just  1.89  seconds  for  a  given 
value  of  weighting  factor  and  Lagrange  multiplier. 
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5 . 5  Analysis  of  Results 

This  study  primarily  examines  a  relatively  small  set  of 
clear  management  plans,  each  somewhat  flexible  to  possible 
changes  in  future  conditions.  In  terms  of  dissolved  oxygen 
content,  the  objective  is  achieved  when  the  deficit  is  kept 
consistently  egual  to  or  close  to  the  acceptable  limit. 
Generally,  this  maximum  allowable  deficit  is  prescribed  by 
some  regulatory  agency.  Hence,  a  system  would  be  called 
good  if  the  critical  deficit  continually  approached  this 
limit,  but  never  exceeded  it.  Whenever  the  water  quality 
goes  below  the  specified  limit,  then  aeration  should  be 
provided  to  keep  the  water  quality  above  the  limit.  Now  the 
objective  is  tc  minimize  the  total  quantity  of  aeration 
supplied  to  maintain  water  quality.  In  this  context  the 
results  indicates  considerable  promise  for  the  use  of 
dynamic  programming  to  achieve  the  specified  water  quality 
standard . 

The  Figure  5.5  and  5.6  are  the  initial  profiles  of  DO 
and  BOD  concentration  for  different  amount  of  pollution  load 
without  artificial  aeration.  This  profiles  are  obtained  by 
solving  equations  4. 1  and  4.2  with  initial  boundary 
conditions.  Figure  5.7  through  5.11  shows  a  typical  optimal 
aeration  capacity  allocation  policy  and  the  corresponding  DO 
sag  profile  for  the  case  €  equal  to  0.1,  weighting  factors 
Wl  =  0.9  and  W2  =  0.  I  and  for  different  quantities  of 


vj  ft 


. 


-  '  - 


.  1  k 


. 


77 


DO  CONCENTRATION  UITHOUT  WASTE  FLOW  AT  REACHES 
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BOD  CONCENTRATION  UITH0UT  UASTE  FLOW  AT  REACHES 


=0.9,  W  =0.1,  MAXIMUM  POLLUTION  (BOD)  =  30  Mg/1. 
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AERATION  CAPACITY  AND  DO  SAG  PROFILE 


=0.1,  MAXIMUM  POLLUTION  (BOD)  =  25  Mg/1 
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pollution  load.  The  solid  line  gives  the  DO  sag  profile 
without  artificial  aeration. 

In  order  to  show  the  effects  of  various  weighting 
factors  on  the  optimal  policy,  the  following  five  sets  of 
weighting  factors  have  been  tested  to  each  of  17  different 
values  of  €. 

(W1,W2)  =  ( • 9,  .  1 )  ,  (.7,. 3),  (.5,. 5),  (.3, .7),  (.1,.9) 

Figure  5.7  shows  a  typical  optimal  aeration  capacity 
allocation  policy  and  the  corresponding  DO  sag  profile  for 
the  case  of  €  equal  to  0.1  and  weighting  factors  HI  =  0.9 
and  W2  =  0.  1.  In  Figure  5.7  the  solid  line  gives  the  DO  sag 
profile  without  using  artificial  aeration.  The  total 
available  aeration  capacities  corresponding  to  different 
values  of  6  and  weighting  factors  are  given  in  TABLE  5.3. 

By  using  the  results  of  TABLE  5.3,  Figure  5.12  was  plotted 
to  illustrate  the  typical  relationship  between  the  total 
available  aeration  capacity  and  different  values  of  6. 

As  mentioned  before,  the  main  concern  of  the  dynamic 
programming  approach,  therefore,  is  to  show  how  artificial 
instream  aeration  can  be  optimized  so  as  to  minimize  the 
relative  cost  of  design  and  operation  for  several  competing 
design  criteria.  With  this  in  view  the  following  results 


were  observed. 


. 


"  L  \ 


.  .  . ,  •  •  • 


. 

•  :  • 

. 


85 


———————  *  Total  Available  Aeration  Capacities 

Corresponding  to  Different  Values  of  6 
and  Weighting  Factors 
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The  amount  of  aeration  required  to  maintain  the 
specified  water  quality  standard  increases  with 
increasinq  waste  discharqe  (Figures  5.7  to  5.11).  For 
a  stream  with  a  BOD  concentration  of  less  than  10 

the  DO  content  of  the  whole  stream  remains  above 
5  mg/1  and  thus  no  aeration  is  required  (Figure  5.10). 
However,  if  the  BOD  level  is  raised  beyond  20  mg/1, 
artificicai  aeration  is  needed  (Figures  5.7,  5.8,  5.9 
and  5.11)  . 

Since  W 1  and  W2  are  relative  weighting  factors,  the 
cost  function  Z  constitute  a  relative  cost  in 
optimization.  The  dimension  of  Z  is  ’dollars'  or 
dimensionless,  depending  on  how  one  interprets  Wl  and 
W2.  In  any  event  the  numeric  values  of  Z  computed  are 
proportional  to  the  true  cost  in  optimization. 

Further,  higher  ratio  of  81/K2  gives  a  better  DO  sag 
profile  in  meeting  the  stream  standard  for  a  given 
value  of  6,  and  one  can  eventually  get  DO  sag  profile 
which  is  completely  below  the  given  stream  standard  by 
choosing  a  very  large  WI/W2  ratio  and  a  proper  e 
value.  DO  sag  profiles  below  the  stream  standard  have 
been  observed  when  6  is  less  than  2  and  W1/W2  is  set 
to  9.  This  shows  that  the  method  can  be  applied  to  a 
pollution  control  problem  in  which  the  stream  standard 
is  required  to  be  satisfied.  In  addition,  the  method 
can  also  be  used  in  a  preliminary  analysis,  since  by 


. 

. 


.1 


.  (II. 


' 

. 


.c 


■ 

’  .1 

. 


88 


varying  the  values  of  6,  different  combinations  of  the 
total  aeration  capacity  used  and  the  resulting  DO 
profile  can  be  obtained. 

3.  The  optimal  total  aeration  capacity  increases  with 
descreasing  6  value.  This  fact  is  clearly  depicted  by 
either  Figure  5.12  or  TABLE  5.3.  As  stated  before  the 
lagrange  multiplier  €  has  the  meaning  of  price,  large 
value  of  6  means  a  higher  cost  for  allocating  an  unit 
aeration  capacity  into  the  system.  The  amount  of 
total  aeration  capacity  will  increase  with  value  of  € 
goes  down  (Figure  5.12).  The  rate  of  increase  will 
decrease  gradually.  The  maximum  capacity  is  attained 
when  6  equals  to  zero.  This  is  the  limiting  case  in 
which  there  is  no  limitation  on  the  total  available 
aeration  capacity. 

4.  The  optimal  total  aeration  capacity  decreases  with  the 
increasing  weighting  factor  K2,  this  implies  that  if 
the  aeration  costs  are  weighted  more  and  more  heavily, 
the  total  available  aeration  capacity  will  decreases. 
This  fact  can  be  observed  from  TABLE  5.3. 

5.  There  exists  a  certain  range  for  the  total  available 
aeration  capacity  with  respect  to  the  change  of  6, 
beyond  this  range,  any  change  in  6,  will  not  affect 
the  amount  of  total  resource.  A  near  optimal  policy 
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can  be  obtained  within  the  gap  by  an  approximation 
technique  (Everett  III  1963).  Figure  5.12  demonstrate 
that  the  range  begins  approximately  when  €  becomes 
greater  than  four. 
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CHAPTER  SIX 


CONCLUSION 


One  of  the  problems  in  river  pollution  control  is  to 
find  an  acceptable  condition-cost  relation  so  that  agreement 
can  be  reached  and  the  necessary  control  established.  In 
this  paper  an  attempt  is  made  to  build  a  model  of  the  system 
which  reveal  certain  salient  features  of  the  system.  The 
main  demand  is  for  interaction  between  water  pollution 
engineer  and  the  computer,  so  that  the  engineer  may  maintain 
overall  control  and  evaluate  the  system  over  intermediate 
results  to  modify  the  parameters  if  he  desires. 

The  method  of  analysis  presented  herein  is  a 
significant  departure  from  the  traditional  methods  for  river 
pollution  control.  In  current  practice,  it  is  customary  to 
go  for  either  a  low  flow  augmentation  dams  or  the  enhanced 
in-plant  treatment.  The  enhanced  in-plant  treatment  is 
widely  used  nowadays;  however,  its  inflexibility  and  its 
relatively  high  cost  do  not  necessarily  make  it  the  best 
method  for  waste  abatement.  The  low-flow  augmentation 
reservoir  is  expensive  and  further,  it  is  net  possible  in 
certain  situations.  Artificial  aeration,  though  still  in 
the  stage  of  development,  offers  many  special  advantages. 

The  use  of  lagrange  multiplier  in  the  solution  of 
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allocating  limited  resources  allows  the  original 
constrainted  optimization  problem  to  be  simplified  to  an 
unconstrained  one.  In  general,  different  values  of  €  lead 
to  different  total  available  aeration  capacities,  hence  one 
can  find  all  possible  available  capacities  and  its 
corresponding  optimal  allocation  policy  by  changing  the 
value  of  €  successively  over  a  wide  range.  Then  a 
parametric  curve  for  the  total  available  aeration  capacity 
with  respect  to  6  can  be  plotted  for  a  given  range  of 
pollution  value.  Once  a  total  aeration  capacity  is  selected 
based  on  the  budget,  one  can  find  the  corresponding  6  value 
using  Figure  5.  12.  By  utilizing  this  value  of  €  in  the 
modified  system,  equations  4. 16  to  4. 19,  one  can  find  the 
optimal  allocation  policy  for  the  given  aeration  capacity  by 
solving  the  discretized  dynamic  programming  algorithm. 

The  use  of  a  dynamic  programming  algorithm  in  this 
nonlinear  optimization  problem  is  very  efficient  in 
obtaining  the  optimal  policy  for  the  given  aeration 
capacity.  Computational  effort  is  the  limiting  factor  in 
the  size  of  the  problem  that  can  be  handled,  tnat  is,  the 
number  of  discretized  states  chosen  depends  on  the  computer 
time  and  storage  available.  Most  of  the  computational 
effort  is  involved  in  the  evaluation  of  objective  function. 
For  a  given  number  of  stages (N)  the  number  of  times  the 
objective  function  to  be  solved  is  directly  proportional  to 
the  number  of  states  in  each  stages.  That  is,  doubling  the 
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number  of  states  will  approximately  double  the  computational 
effort.  In  this  example  there  were  24  stages  with  I  I 
discretized  states  at  each  stage  and  the  dynamic  programming 
algorithm  find  the  optimum  value  out  of  all  II24  possible 
aeration  policies  in  just  1.89  seconds  for  a  given  value  of 
weighting  factors  and  lagrange  multiplier. 

The  existance  of  optimal  value  for  on  objective 
function  of  the  discretized  type  variable  is  assumed  by  the 
algorithm  itself,  since  all  possible  ways  have  been  evaluted 
and  compared  by  the  algorithm.  Philosophically  speaking,  it 
is  generally  impracticable  to  implement  the  entire  optimal 
solution  in  practice;  in  fact  it  is  seldom  necessary  to  do 
this.  This  optimal  control  solution  can  serve  a  very  useful 
purpose,  however,  in  establishing  limits  of  performance  and 
relative  cost.  It  is  perhaps  even  more  useful  in  comparing 
the  optimal  design  policy  with  several  competing  suboptimal 
designs. 

The  model  described  in  this  work  is  limited  to  those 
needed  to  demonstrate  the  general  utility  of  dynamic 
programming  approach  to  river  pollution  control.  Due  to  the 
comprehensive  nature  of  the  model  and  its  flexible 
procedures,  many  other  input-output  studies  are  possible. 

The  model  is  inherently  quite  flexible,  and  any  or  all  of 
the  relationships  may  be  modified  without  changing  the 
dynamic  programming  methodology. 
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It  is  emphasized  that  the  choosen  relationships  are  not 
necessarily  the  best  for  every  day  and  every  season.  But 
the  concept  of  dynamic  programming  approach  for  the  control 
of  river  pollution  is  most  significant  in  this  paper. 
Methodology  such  as  the  one  the  author  has  developed  is 
useful  for  the  present,  only  at  an  elementary  stage.  Much 
research  is  needed  before  its  use  could  become  practicable. 
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APPENDIX-A 


OSEEVS  MANUAL 

The  material  contained  in  this  appendix  is  intended  for 
those  making  the  detailed  preparation  for  the  use  of 
CHKDATA,  RIPSIM  (River  Pollution  Simulation  Model)  and 
OPTMOD  (OPTmization  MODel) .  This  appendix  has  been  prepared 
in  sufficient  detail  to  serve  as  a  manual  for  guidance  in 
the  preparation  of  program  inputs  and  control  statements. 

The  information  is  separated  into  program  components  and 
each  component  is  described  in  the  following  pages. 
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FIGURE  A .  1 
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A 1 .  1  .C  HKD  AT  A  -  Stream  Flow  Data  Edit  Program 

£il££^se:  The  CHKDATA  program  is  designed  to  read  daily  raw 
streamflow  data  of  the  type  available  from  the  Federal  Water 
Resources  Branch,  Edmonton,  Alberta,  for  their  stream  gaging 
stations  and  to  prepare  these  data  for  use  in  synthetic 
streamflow  data  generator.  The  program  reads  the  raw  data 
from  magnetic  tape,  checking  as  it  does  so  to  ensure  that 
the  data  being  read  are  in  the  proper  sequence,  for  the 
proper  station  and  for  the  years  of  record  desired. 

When  all  desired  raw  data  for  one  month  have  been  read  in 
the  proper  order,  the  program  searches  for  missing  data.  If 
more  than  one  month  of  consecutive  data  are  missing,  the 
program  calls  exit  and  the  operator  must  reschedule  the  data 
sequence  to  remove  that  period  from  use.  If  fewer  than  one 
month  of  data  are  missing,  the  program  fills  the  missing 
data  so  that  the  output  is  a  complete  set  of  daily 
streamflow  records. 

Missing  data  are  filled  by  noting  the  day  of  the  month  and 
year  of  the  missing  data  item  and  its  station  number.  Then 
a  mean  of  all  other  days  of  that  month  data  item 
corresponding  to  the  one  which  is  missing  is  computed.  The 
mean  of  mean  and  standard  deviation  for  that  month  are 
computed  from  these  data  items  and  then  the  missing  item  is 
computed  according  to  the  formula: 
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Q  (i)  =  mean  (i)  +  std(i)  *  r  (i) 
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Where 

Q(i) 
mean  (i) 

std  (i) 


r(i) 


the  computed  missing  data  item, 
the  mean  of  all  data  available 
corresponding  to  that  month, 
the  standard  deviation  of  all  data 
available  corresponding  to  the  month  and 
station  of  the  missing  data  item, 
a  standard  normal  random  deviate. 


Having  read,  edited  and  filled  the  missing  data,  the  program 
then  computes  and  outputs  the  mean  flow  and  standard 
deviation  for  that  month.  There  are  12  means  and  12 
standard  deviations  for  each  station.  These  means  and 
standard  deviations  are  used  in  the  synthetic  data 


generator . 
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A 1 . 2  Program  Components 

The  CHKDATA  program  consists  of  the  subroutines,  listed  with 
their  lengths  in  bytes,  as  follows: 


CHKDATA  MAIN -  105  10  A VM - 706 

INPOT  1288  RAN - 832 

INCARD  1640  RANDO  —  448 

FILL  1332 


CHKDATA  MAIN 

This  program  component  reads  the  program  contolling 
information,  coordinates  the  work  of  the  other  subroutines, 
makes  certain  checks  and  outputs  the  edited  data.  The 
program  controlling  information  is  supplied  on  two  cards 
(explained  in  Program  Input  and  Output  section) .  The 
controlling  information  establishes  the  waste  disposal 
points,  the  year  in  which  data  for  each  point  begin  and  end. 

After  identifying  the  points,  the  subroutines  INPUT  and 
INCARD  are  called.  These  subroutines  read  the  data  one 
month  at  a  time,  perform  certain  checks  described  below,  and 
return  for  one  month  to  CHKDATA  MAIN.  Checks  for  the 
station  number  and  beginning  and  ending  years  then  are  made 
on  the  month’s  data  just  read.  Any  deviation  from  proper 
order  is  printed. 
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When  the  reading  of  data  is  completed,  the  subroutine  FILL 
is  called  to  fill  in  missing  data.  FILL  is  described  below. 
Finally,  depending  upon  the  controlling  information 
supplied,  CHKDATA  MAIN  outputs  the  edited  streamflow  data. 
The  data  are  supplied  in  printed  form,  or  on  magnetic  tape. 
Daily  flows  are  outputed  without  further  change  (if 
required) .  Monthly  flows  are  the  average  of  the  daily  flows 
for  one  month. 

subroutine  INPUT 

This  subroutine  transmits  the  call  to  read  data  from 
CHKDATA  MAIN  to  subroutine  ICARD  and  the  data  read  back  to 
CHKDATA  MAIN.  Before  the  data  are  sent  to  CHKDATA  MAIN, 
this  subroutine  checks  to  determine  that  the  data  read  are 
of  proper  station,  month  and  year.  The  data  are  tramsmitted 
in  a  one  dimensional  array,  one  month’s  data  at  a  time. 

Subroutine  ICARD 

ICARD  actually  reads  the  supplied  raw  historical  data. 
The  subroutine  checks  to  determine  if  the  data  are  in  the 
proper  monthly  sequence.  Data  for  one  month  at  a  time  are 
read,  checked,  and  transmitted  to  subroutine  INPUT. 

Subroutine  FILL 

After  the  data  are  all  read  and  checked,  subroutine  fill 
is  called  to  determine  if  there  are  any  missing  data  points. 
If  thirty  or  more  consecutive  daily  data  points  are  missing 
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for  any  station,  exit  is  called  and  the  operator  must  make 
an  adjustment  in  the  data  years  used.  If  scattered  data 
points  are  missing,  the  subroutine  fills  the  missing  points 
one  at  a  time,  as  described  above. 

Subroutine  RAN  and  RANDU 

These  subroutines  generate  standard  normal  deviates 
(mean=0,  and  variance=l)  for  use  in  subroutine  FILL. 

Subroutine  AVM 

Subroutine  AVM  compute  the  monthly  average  flows  from 
the  edited  and  filled  daily  data  and  outputs  on  a  magnetic 
tape,  disk  or  cards,  depending  on  the  control  number  entered 
for  the  variable  ITAPS. 
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A  1 . 3  Program  Input  Form  at 


Input  format  and  the  meaning  of  certain  important 
variables  are  as  follows: 


CHKDATA  MAIN 
Card  #1  FORMAT  (915) 

IYR 1  (I) 
IYR2  (I) 
NNSTA 
ITAPE 

ISTART 

ITIME 

IPRINT 

NTAPE 


=  The  starting  year  for  station  (I) . 
=  The  ending  year  for  station  (I) . 

=  The  number  of  disposal  points. 

=  4  for  data  input  on  tape. 

=  5  for  data  input  on  cards. 

=  A  starting  random  number. 

=  1  for  daily  flow  output. 

=  2  for  average  monthly  flow  output. 
=  0  if  output  is  not  to  be  printed. 

=  I  if  output  is  to  be  printed. 

=  0  if  output  is  not  to  be  taped 
=  1  if  output  is  to  be  taped. 


Card  #  2  FORMAT(10I8)  (ISTA(I)  ,  1=  1 ,  NNSTA) 


Indicates  the  station  number  for  which  data  are  to  be 
read.  These  station  numbers  must  correspond  to  the 
reach  numbers.  The  number  of  stations,  NNSTA,  may  be 
fewer  than  the  number  of  reaches  but  the  order  in 
which  stations  are  read  must  be  the  same  as  the  reach 


numbers. 
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Subroutine  ICARD 

A  set  of  input  streamflow  data  cards  are  required  to  input 
all  stream  flow  data.  The  data  format  is: 

FORMAT  (18,14,12, 8F6.2) 


where. 

18 

=  the 

station 

number. 

14 

=  the 

calender 

year. 

12 

-  an 

Index,  1 

,  2,  3,  or  4,  which 

identifies  the  number  of  the  card 
in  the  month. 

8F6 . 2  =  The  daily  streamflow  values  for  8 

consecutive  days.  The  first  card 
contains  flow  data  for  the  first  8 
days,  the  second  for  the  9th 
through  the  16th  day,  the  third  for 
the  17th  through  the  24th  day  and 
the  fourth  for  the  25th  day  through 
the  last  day  of  the  month,  28,29,30 
or  31,  as  appropriate. 
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A  I  .  4  Dictionary  of  Variables. 

Following  is  a  list  of  variables  used  in  C HKD AT A  and  their 


meanings 

in  brief: 

AV2  (I, J) 

Average  flow  daring  Ith  month  of  the  Jth  year 

FN 

Number  of  data  items. 

ICARD(I) 

Card  sequence  number. 

ICOUNT 

Convenience  Index. 

IDAY 

Day  counting  variable 

IISTA 

Convenience  station  number,  for  checking. 

IIYRI 

Convenience  starting  year,  for  checking. 

IIYR2 

convenience  final  year,  for  checking. 

IPUNCH 

Control  variable-output. 

ISTART 

Starting  random  number. 

1ST  A  (I) 

Identifying  number,  Ith  station. 

ITAPE 

Control  variable  for  output, 

=4  data  on  tape 
=5  data  on  cards 
=6  data  on  disk  file 

IT  I  ME 

Control  variable-input, 

=1  for  daily  flows. 

=2  for  monthly  average  flows. 

I MO NTH 

Month  counting  variable. 

IYEAR 

Year  counting  variable. 

IYR 1  (I) 

Starting  year  for  data,  Ith  station. 

IYR2  (I) 

Final  year  for  data,  Ith  staion. 

JPREV 

Convenience  station  number,  for  checking 

JYR 

Convenience  number  of  years  for  checking. 

KMO  (I) 

Identifying  number,  Ith  month. 

. 

. 
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KSTA  (I) 

Equals  ISTA  (I)  . 

KYF.  (I) 

Equals  IYR  (I)  . 

LHC 

Number  of  months  of  data,  computed. 

MONTH 

The  current  month. 

ND  (I) 

Number  of  days  in  month  I 

NEVEN 

Convenience  number,  random  number  generator. 

NMO 

Number  of  months  of  data  counted. 

NNSTA 

The  number  of  stations. 

NNYR 

Number  of  years  of  data. 

HODD 

Convenience  number,  random  number  generator. 

NSTA 

Current  station  number. 

NTAPE 

Control  variable-output, 

=0  for  no  data  on  tape. 

=1  for  data  output  on  tape. 

Q(I,J,K) 

Gage  flow  Ith  month,  Jth  day,  Kth  year. 

QMEAN 

Mean  flow. 

QNEW 

Computed  flow  to  fill  missing  data  point. 

QSTD 

Standard  deviation  of  flow. 

R 1 ,R2, YFL 

Convenience  variables,  random  number 
generator . 

RN  (I) 

Random  number,  Ith  time  frame. 

SUM 

Sum  of  second  term. 

S(I,J) 

Flow  for  Jth  day  of  Ith  month. 

Z(I) 

Temporary  storage,  flow  on  Ith  day  in  current 

month. 

Z(L) 

Temporary  linear  storage  for  one  month* s  flow 

. 
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A2. 1  RIPSIM  (Elver  Pollution  Simulation  Model) 


The  simulation  program,  called  RIPSItM,  is  made  up  of  a 
controlling  program,  RIFSIM-MAIN  and  several  subroutines 
each  of  which  contributes  to  the  overall  program. 


A 2 . 2  Program  Components . 

A  list  of  the  subroutines  and 
given  below. 


RIPSIM  main 

906 

SIM 

1544 

UPGAGE 

1  162 

OPTMOD 

15004 

RAN 

560 

RRN 

396 

STD 

490 

RIPSIH  MAIN. 


their  length,  in  bytes  are 


1 WASTE 

522 

TGEN 

2772 

TRAN 

2852 

IREACH 

4  12 

COM  PUT 

8580 

GTRAN 

536 

GFLOK 

456 

Purpose:  The  simulation  program  is  controlled  by  RIPSIM 
main.  This  routine  sets  up  ten  common  blocks  needed  to 
transfer  variable  values  from  subroutine  to  subroutine, 
calls  subroutine  in  the  required  order  and  reads  the  station 
numbers  to  be  used  in  simulation.  Subroutines  TGEN  and  TRAN 
are  called  to  set  up  the  reach  indices  and  computation 
sequence  and  to  compute  the  reach  flow  data.  The  subroutine 
SIM  is  called  to  carry  out  the  simulation  process. 
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Subroutine  SIM 


Subroutine  SIM  reads  in  constants  and  data  needed  to  compute 
the  temperature  and  waste  parameter  values  for  each  reach, 
computes  the  deterministic  component  of  the  temperature 
equation,  and  through  two  independent  "do  loops"  to  perform 
simulation  by  calling  various  subroutines.  The  first  of  the 
"do  loop”  covers  from  the  first  reach  to  N,  the  number  of 
reaches.  This  "do  loop”  also  calls  RAN  which  generates  the 
random  deviates  for  the  temperature  equation. 

The  second  "do  loop"  computes  DO  deficit  and  checks  for 
water  quality  violation.  Inside  this  "do  loop"  subroutine 
QTRAN  and  OPTMOD  are  called,  successively,  to  compute  flow 
and  the  water  quality  values  at  each  reach.  If  there  is  a 
violation  in  water  quality,  the  subroutine  OPTMOD  is  called 
to  compute  the  amount  of  aeration  required  to  avoid 
violation.  These  operations  are  described  below  in  more 
detail  as  individual  subroutines.  The  values  of  the 
simulation  results  are  written  out  within  the  "do  loop". 

When  the  simulation  is  complete,  the  control  will  be 
transfered  to  RIPSIM  main  for  termination. 

Subroutine  RAN  and  E ANDU 

Subroutines  RAN  and  RANDU  are  the  same  random  number 
generating  subroutines  used  in  program  CHKDATA  which  have 
been  described  previously. 
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Subroutine  GFLOW 


Subroutine  GFLOW  computes  total  flow  at  the  beginning  of 
reach-  GFLOW  is  called  in  the  first  "do  loop"  to  compute 
total  flow  in  the  stream  at  the  beginning  for  each  waste 
disposal  stations  being  used  in  the  simulation.  The 
computed  data  are  transferred  to  COM  MON/FLO  W7  for  use  by 
other  subroutines. 

Subroutine  QTRAN. 

QTRAN  is  a  short  subroutine  which  receives  computed  flow 
data  through  C0MM0N/FL0W7  and,  convert  the  computed  flow 
data  to  stream  flow  data  by  linear  transformation.  The 
streamflow  data  are  designated  QNAT (I) ,  to  indicate  the  flow 
at  the  upper  end  of  reach  I.  These  values  of  QNAT  (I)  are 
also  entered  in  C0MM0N/FL0K7  for  use  by  other  subroutines 
linked  through  this  common  block. 

Subroutine  TWASTE 

TWASTE  checks  whether  there  is  a  set  of  waste  load  data 
corresponding  to  a  reach  number  where  a  waste  discharge  is 
scheduled.  If  the  data  are  not  available  for  the  reach 
TWASTE  calls  RANDU  to  generate  waste  load  for  that 
particular  reach. 
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Subroutine  OPTMOD 


After  regulated  flows  in  each  reach  are  computed,  subroutine 
OPTMOD  is  called  to  compute  the  water  quality  values  for 
each  reach  beginning  from  the  first  reach.  First,  OPTMOD 
completes  the  computation  of  temperature  for  the  time-frame. 
It  should  be  remembered  that  OPTMOD  is  called  in  "do  loop" 
which  is  in  turn  a  "do  loop".  OPTMOD  sets  up  a  third  "do 
loop"  which  cycles  over  the  number  of  reaches,  NR.  Thus, 
for  each  reach  ,  OPTMOD  compute  the  quality  values  for  each 
reach  before  it  returns  to  its  calling  subroutine,  SIM. 

The  computations  are  made,  starting  at  the  downstream 
reaches,  according  to  the  seguance  of  computation  set  by 
subroutine  TGEN. 

For  each  reach,  OPTMOD  initilizes  variables  and  proceeds  to 
make  the  computations  necessary  to  evaluate  the  incoming  BOD 
and  DO  concentrations,  K1  (deoxygenation  constant),  K2 
(reaeration  constant)  and  corrects  K1  and  K2  for 
temperature,  computes  the  time  of  flow  in  each  reach  and, 
finally,  the  BOD  and  DO  concentrations  in  water  leaving  the 
reach.  OPTMOD  then  checks  to  determine  if  the  DO  deficit 
has  reached  a  maximum  within  the  reach,  in  which  case  the 
minimum  DO  concentration  will  have  occured  within  the  reach. 
This  is  done  by  computing  TCRIT,  the  critical  time  of  flow, 

and  comparing  it  to  the  time  of  flow  in  the  reach.  If 

DO  concentration  has  occurred  in  the 


TCRIT  <  TIME,  a  minimum 
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reach.  TCRIT  is  then  substituted  for  TIME  and  the  value  of 


DO-OUT  is  computed  and  used  in  turn  to  compute  the  minimum 
DO  concentration,  XMINDO,  in  the  current  reach.  If  there  is 
any  violation,  then  optimization  program,  OPTMOD,  computes 
artificial  aeration.  OPTMOD  then  writes  out  the  water 
quality  and  related  values  in  an  array  and  then  returns  to 
its  calling  routine. 
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A  1 . 3  Program  Input 


1  1  7 


Data  are  read  into  the  simulation  program,  RIPSIM,  through 
RIPSIM  main,  SIM,  GF10W,  and  TGEN  only. 

RIP SIM  MAIN 
Card  #1  FORMAT  (215) 

=  the  number  of  times  of  simulation  to  be 
carried  out. 

NGT  =  the  number  of  stations  for  which  data 

are  used. 


Card  #2  FORMAT  (218) 

IGT(I)  =  the  waste  disposal  station  numbers  used. 

I  =1,....,  NGT.  Maximum  of  twenty  four. 

Subroutine  SIM 


Card  #1  FORMAT (3F8 . 0 ,  15,  F8.0) 


D 

C 

TMEAN 

ISTART 

SIGMAT 


=  constant  in  temperature  eguation. 

=  lag  constant  in  temperature  eguation. 

=  mean  annual  temperature. 

=  initial  number  for  random  number  generator, 
=  standard  deviation  of  the  temperature  data, 


Card  #2  FORMAT  (2F 1 0 . 0 ,  II) 


RLNTH (I) 
XK  1  20W  (I) 


=  the  length  of  reach  (i) ,  in  days. 

=  the  deoxygenation  constant  K  I  at  20°C 
for  the  waste  which  is  introduced 
into  reach  I. 


IWASTE  (I) 


=1  if  there  is  a  waste  load  introduced 


in  reach  I. 
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-0  if  no  waste  load  introduced  in  reach  I. 
Card  #3  FORMAT  (3F  I  0) 

QWASTE(I,J)  =  the  rate  of  discharge  of  waste  in  Cfs  for 

reach  I  and  month  J. 

BODWST  (I, J)  =  the  BOD  concentration  in  the  waste  being 

discharged  into  reach  I  during  month  J. 
DOWST(I,J)  =  the  DO  concentration  in  the  waste  being 

discharged  into  reach  I  during  month  J. 

Subroutine  GFLOW 

The  computed  data  are  read  by  GFLOF.  The  format  must  be 
adjusted  to  the  format  of  the  data  to  be  read.  Normally, 
the  data  will  be  on  magnetic  tape  .  The  program  as  set 
forth  herein  reads  the  data  from  magnetic  tape  as  QG(I,J), 
(1=1,  NG)  ,  (J=  I,  12).  It  reads  one  year  of  data  for  each 

station. 

Subroutine  TGEN 
Card  #1  FORMAT  (15,  FI0.0) 

NR  =  number  of  reaches. 

NG  =  number  of  waste  disposal  stations. 

Card  #2  One  card  for  each  reach  (215) . 

NOR  (I) 


FL  (I) 


=  the  reach  number,  reach  I. 
=  length  of  reach  I. 
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A. 2 . 4  Dictionary  of  Variables 


A 

Intermediate  constant. 

AVW (I) 

Average  waste  load,  Ith  reach 

(used  for  periodic  load  function) . 

AW  (I) 

Amplitude  of  waste  load  (used  for  periodic 

load  function. 

E 

Intermediate  constant. 

BODIN  (I) 

BOD  concentration  at  the  end  of  reach  I. 

BODCUT  (I) 

BOD  concentration,  leaving  reach  I. 

BODWST  ( I) 

BOD  concentration  in  waste  for  reach  I. 

C 

constant,  temperature  eguation. 

CN 

correction  value,  flow  in  upstream  reach. 

CNI 

Intermediate  variable. 

CQREG(I) 

corrected  regulated  flow,  reach  I. 

2(1,  J) 

A  state  variable  at  stage  I, 

state  position  J. 

DEFIN  (I) 

DO  deficit  concentrati on ,  upstream  end 

of  reach  I. 

DEFOUT  (I) 

DO  deficit  concentration,  leaving  reach  I. 

DOS 

Dissolved  oxygen  saturation  concentration. 

DOWST  (I,  J) 

DO  deficit  concentration  in  waste  load, 

reach  I,  month  J. 

DQ(I) 

Amount  of  flow  regulation  reach  I. 

F(I,J,K) 

The  up-to-date  total  cost  for  k 

keeping  at  the  present  state  from  stage  1 

stage  I. 

FL  (I) 

Length  of  reach  I. 

. 
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IWASTE  (I) 

Equals  1  for  a  waste  load  in  reach  I. 

Equals  0  for  no  waste  load  in  reach  I. 

JE(I) 

Downstream  reach  index. 

NOE  (I) 

Number  of  reach  I. 

NE 

Number  of  reaches. 

NWASTE 

Number  of  waste  load  stations. 

NYE 

Number  of  years. 

PMW (I) 

Peak  waste  load  in  reach  I  (used  for  periodic 

load  function) 

QG  (I,  J) 

Weekly  generated  flow,  station  I,  month  J. 

QNAT  (I) 

Natural  (unregulated)  flow,  reach  I. 

QEEG  (I  ,  J) 

Eegulated  flow,  reach  I,  month  J. 

QSUM 

Sum  of  incoming  flows,  correc ted-waste  loaded 

reaches  only. 

QSUM2 

Sum  of  incoming  flows,  total. 

CWASTE  (I,  J) 

Sate  of  waste  discharge,  reach  I,  month  J. 

E(I) 

Bandom  number,  normal  distribution, 

time  frame  I. 

BOON  (I) 

Deoxygenenation  error  constant,  reach  I. 

ELNTH  (I) 

Length  of  reach  I. 

EEN  (I) 

Temporary  variable,  random  number 

generator,  time  frame  I. 

EEV 

Intermediate  variable. 

BTEMP 

Intermediate  variable. 

SDIF 

Intermediate  varaiable. 

SIGMAT 

Standard  deviation,  temperature  data. 

SPED  1 

Summing  variable  -  K1  computation. 
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SPRD2 

Summing  variable  -  BOD  computation. 

SPRD3 

Summing  variable  -  DO  deficit  computation. 

T  (L) 

Temperature  for  reach  L. 

TCRIT 

Critical  time  of  flow,  time  to  critical  DO 

condition . 

TIME 

Time  of  flow  in  reach. 

THEAN 

Mean  temperature. 

TT 

Temperature  at  current  time  frame. 

V(I,J,K) 

A  decision  variable  at  state  I. 

It  represents  the  aeration  capacity  needed 

at  point  1-1  when  the  state  variable  at 

stage  I-l  is  D(I-I,J)  and  at  stage  I 

is  D  (I ,  K)  . 

XK120  (I) 

Deoxygenation  constant,  reach  I, 

at  2 0°C . 

XK2 

Reoxygenation  constant. 

XK220 

Reoxygenation  constant  at  20°C. 

XKK  (I) 

Dissolved  oxygen  concentration,  violation 

condition . 

XMINDO 

Minimum  dissolved  oxygen  concentration. 
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APPENDIX-B 


THEORY  OF  OXYGEN  SAG  CURVE 


One  of  the  common  and  important  indicators  of  water 
quality  is  Dissolved  Oxygen  (DO) .  The  oxygen  deficiency  in 
a  stream  at  any  point  is  the  resultant  of  two  opposing 
processes,  each  governed  by  a  different  set  of  conditions  - 
one  consuming  oxygen  from  the  stream  and  the  other 
replenishing  the  oxygen  source,  both  processes  proceeding  as 
a  time  function  (figure  B.  1)  . 

The  principal  wastes  which  affect  DO  are  those  which 
are  characterized  by  their  Biochemical  Oxygen  Demand  (BOD) . 
BOD  is  generally  defined  as  the  amount  of  oxygen  (DO) 
required  by  bacteria  while  stabilizing  decomposable  organic 
matter.  The  stabilization  of  waste  matter  using  dissolved 
oxygen  is  called  purification  of  water. 

River  pollution  control  requires  careful  consideration 
of  the  relationship  between  BOD  and  DO.  One  of  the  primary 
considerations  in  determinig  the  amount  of  waste  (BOD)  which 
can  be  released  to  a  stream  is  the  DO  content  of  the  stream 
at  the  point  of  discharge.  Down  stream  from  this  point  the 
DO  content  is  reduced  by  the  oxidation  of  organic  material. 
At  some  point  downstream  the  rate  of  oxygen  recovery, 
(reaeration)  will  exceed  that  of  oxygen  reducrion 
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(depleation)  and  from  that  point  onwards  the  DO  content  will 
gradually  increase.  The  phenomena  of  depletion  and  recovery 
of  oxygen  level  in  a  stream  can  be  discovered  by  the 
Streeter-Phelps  equations. 


dD 

-  =  k  IL  -  k2D 

Dt 

This  basic  differential  equation  states  that  the  net  rate  of 
charge  in  the  DO  deficit  (dD/dt)  is  equal  to  the  sum  of  (I) 
the  rate  of  oxygen  utilized  by  BOD  in  the  absence  of 
reaeration  (K1L)  and  (2)  the  rate  of  oxygen  absorption  by 
reaeration  in  the  absence  of  BOD  (-K1D). 

The  critical  oxygen  deficit,  Dc,  is  defined  as  the 
dissolved  oxygen  deficit  at  the  time  Tc,  the  oxygen  content 
starts  to  increase,  as  shown  in  figure  B.l.  Frequently  Dc 
is  prescribed  by  regulatory  agencies  to  be  less  than  some 
maximum  allowable  level,  Dl.  If  at  any  time  the  value  of 
Dc  <  Dl,  then  the  river  is  unpolluted.  If  Dc  >  Dl,  then  the 
river  is  polluted.  To  avoid  this  being  happen,  more  air 
should  be  pumped  into  the  river  through  aerators.  This 
phenomena  of  artificially  increasing  the  content  of  DO  is 
called  artificial  aeration.  How  much  air  should  be  pumped 
into  the  the  river  system  to  avoid  pollution  is  being  the 
subject  of  this  paper. 
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SATURATION  LEVEL  OF  DISSOLVED  OXYGEN 
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APPENDIX-C 


NOTATIONS  AND  ABBREVIATIONS  USED  IN  THIS  THESIS . 


Ai 

,  th 

Aerator  located  at  the  top  of  i  reach. 

BOD 

Biochemical  xygen  Demend. 

BODin  = 

Quantity  of  waste  already  in  the  system  which  is 
input  to  the  next  reach. 

B0D  dis  = 

Quantity  of  waste  discharged  at  the  top  of  reach!. 

Cfs 

Cubic  feet  per  second. 

o 

o 

II 

Initial  level  of  Dissolved  Oxygen  (DO)  content  at 
the  top  of  a  reach. 

C1 

minimum  level  of  Dissolved  Oxygen,  specified  by 
the  stream  standard  law. 

cs  = 

Saturation  level  of  DO,  a  function  of  temperature. 

n 

rt 

II 

level  of  Do  at  any  point  t  days  from  the  top  of  a 
reach. 

Cc 

Critical  dissolved  oxygen  deficit. 

DO 

Dissolved  oxygen  content  in  water  (mg/1) . 

II 

o 

Q 

Initial  Do  deficit  at  the  top  of  a  reach  = 

<cs  *  C0) 

El 

Maximum  DO  deficit  specified  by  the  regulatory 
agency. 

DO  deficit  at  any  point  t  days  from  the  top  of  a 
reach  (Cg  -  Ct)  . 

L0 

initial  BOD  content  at  the  top  of  a  reach  =  (EODin 
+  BCDdis  )  * 

EOE  content  at  any  point  t  days  from  the  top  of  a 
reach. 

II 

H 

Deoxygenation  coefficient  (rate  of  consumption  of 
oxygen) . 
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Natural  reaeration  coefficient  (rate  o£  recovery 
of  oxygen. 

miligrams  per  litre  (unit  of  measurement) . 

Numfcer  of  reaches  (in  this  case  N=24) . 

Distance  from  the  top  of  a  reach  expressed  in  time 
units  (days)  =  (distance/velocity) . 

Increment  in  the  level  of  DO  due  to  artificial 
aeration  augmented  at  the  top  of  a  reach. 

Length  of  reach  i  =  t  -  t  .  ,) 

i  l-l 

The  value  of  ohjactive  function. 

Lagrange  multiplier. 
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